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CORRECTION FACTOR TECHNIQUES FOR 
IMPROVING AERODYNAMIC PREDICTION METHODS* ** 

By Joseph P. Giesing, Terez P. Kalman 

** 

and William P. Rodden 
Douglas Aircraft Company 

SUMMARY 

This report describes a method for correcting lifting surface theory so 
that it reflects known experimental data. Specifically the theoretical press- 
ure distribution is modified such that imposed constraints are satisfied 
(e.g., lift, moment, etc.) while minimizing the change to the theoretical 
pressure distribution. It is assumed that a finite element or discretized 
lifting surface method is used, such as either the Doublet or Vortex Lattice 
Methods. 

There are several ways in which the theoretical pressures are modified. 
One is a direct application of a set of correction factors to the pressures. 
This is accomplished by premultiplying the pressures with a diagonal matrix 
of correction factors. A second approach to correcting the theory is to 
modify the downwash. This modification can be accomplished by either multi- 
plying the downwash by a diagonal matrix of correction factors or by adding 
an incremental downwash (which is proportioned to the pressure) to the theo- 
retical downwash. In any case the correction factors are adjusted so that 
the imposed experimental constraints are satisfied by the corrected pressure 
distribution while the changes in the pressure distribution are minimized. 

There are several features that have been built into the basic method 
and these include: (1) the ability to consider together experimental data 


* The authors wish to acknowledge Dr. Edward Albano for interesting 
discussions of alternate formats (non-diagonal) of correction matrices and 
their potential derivations. 

** Consulting Engineer 
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from more than one mode (e.g., control surface rotation, pitch, camber, etc.), 
(2) the ability to limit the excursions of the correction factors (i.e., 
establish minimum and maximum values for them) and (3) the ability to use 
correction factor mode shapes (i.e., construct correction factors from known 
distributions or functions). 

The methods developed have been implemented on the computer and many 
correlations and calculations made. Specifically cases involving all three 
Mach Number ranges are considered. For instance in the subsonic speed range 
a swept wing with an oscillating partial span flap and a swept wing with a 
leading edge droop are discussed. In the transonic speed range a two-dimen- 
sional symmetric airfoil with an oscillating flap is treated in detail. An 
arrow wing with and without camber is used in the supersonic analysis. 

The computer program used to generate the correction factors for these 
cases is also fully described in this report and test cases are provided. 

Finally, a new, simple method for accounting for transonic effects in 
the lifting surface theory is described and correlated for the two-dimensional 
case. Basically, a transformed distance between the sending element and re- 
ceiving point is employed. The transformation depends on the time delay en- 
countered by a signal traveling from the sending point to the receiving point. 
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INTRODUCTION 


Wind tunnel data have provided the basis for semi -empirical methods of 
aeroelastic analysis for many years, whether in the estimation of stability 
and control characteristics, the calculation of structural loads, or in flutter 
analysis by modified strip methods. These semi -empirical methods have been 
tailored to aerodynamic lifting-line theory or to strip theory and not to the 
more general (and more accurate) lifting-surface methods. The use of a 
diagonal correction matrix to be applied as a premultiplying factor to matrices 
of aerodynamic influence coefficients obtained from lifting surface theory has 
been considered by a number of authors. A premultiplier may be regarded as a 
correction to the pressure distribution; as an alternative, a postmultiplier 
would be regarded as a correction to the downwash to account for thickness 
effects and for camber induced by boundary-layer displacement effects. Rodden 
and Revel 1 (refs. 1-4) considered a real correction matrix derived from 
static wind tunnel measurements and theoretical load predictions. Bergh and 
Zwaan (refs. 5 and 6) investigated a complex correction matrix derived from 
oscillatory wind-tunnel pressure measurements and theoretical predictions. 

These authors assumed measurements were available only for a single mode, a 
steady angle of attack or an oscillatory pitching (or yawing) mode. 

Current interest in using actively controlled aerodynamic surfaces to 
minimize aeroelastic response requires an improvement in accuracy in predict- 
ing unsteady aerodynamic characteristics of lifting surfaces equipped with 
control surfaces. The correction matrix provides one means of improving the 
accuracy but it requires experimental data on control surface characteristics 

in addition to the angle of attack characteristics of the surface. Hence, an 
extension of references 1 - 4 is necessary to obtain the correction matrix 
for more than one aerodynamic mode. Furthermore, the discrepancies between 
theory and experiment in predicting trailing-edge control surface loads are 
most likely caused by boundary-layer displacement effects on the effective 
downwash. Hence, another extension is necessary to obtain a postmultiplying 
correction matrix. These two extensions are considered in the present 


3 



development. The diagonal format has been retained and complex pre- and . 
postmulti plying correction matrices have been derived which satisfy the con- 
straints of matching experimental data from multiple aerodynamic downwash modes. 

The use of correction matrices is in the time-honored engineering tradition 
of empirical correction factors. It retains the generality of the theory 
while approaching the limiting values of the test results.. Such a posteriori, 
adjustment obviously cannot be regarded as addressing any of the fundamental 
causes of the discrepancies. Other possibilities exist whereby empirical 
corrections can be introduced directly into the theoretical solution. Ashley 
(ref. 7) has discussed two such "irrational correction methods". The first 
of these is of interest here and concerns the calculation of the downwash 
boundary condition and the pressure distribution by "local linearization" 
in terms of the local velocity V L rather than the free stream velocity U m . 

The dimensionless downwash then becomes 

l !x,y ’ 0,t) = r[f + v L (x - y) t] 

where h(x,y,t) is the deflection of the mean surface. Applications of this 
"local linearization" to the downwash boundary condition (but not to the 
kernel function nor pressure coefficient) have been made for control surfaces 
by Ashley and Rowe (ref. 8) and by Rowe, Winther, and Redman (ref. 9), and 
improved correlations have been obtained. Tijdeman and Zwaan (ref. 10) have 
also employed the local linearization of the downwash boundary condition but 
have suggested another modification for use in the Doublet-Lattice Method 
for high subsonic flows, viz., that the free stream Mach number be replaced 
by a mean Mach number for each panel (lifting element or box). The down- 
wash induced at box i by the lifting pressure on box j then becomes 


* The format of a full matrix was briefly investigated but it was found 
to destroy the distributional character of the theoretical aerodynamic in- 
fluence coefficients, and was not considered further. 
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D ij (M jl’ k r )AC pj 


U 

where D.. • (M^, k f ) is the downwash influence coefficient between the jth 
lifting element and the ith downwash collocation point and its functional 
dependence on M^ and the reduced frequency k p is indicated. Tijdeman and 
Zwaan denote the freestream Mach number by M^ and the local Mach number at the 
surface of box i by M i2> so that the locally linearized downwash for harmonic 

motion becomes 


The values of M ] and M 2 for a certain box are not equal, in general, because 
M-j has to reflect the influence of the Mach number distribution normal to the 
surface ranging from M 2 at the surface to the freestream Mach number far 
away from the wing. Preliminary results from NLR calculations have shown 
that M-| can be chosen simply to be the average value of M 2 and M^. 
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LIST OF SYMBOLS 


Matrix that gives pressures in terms of downwash. 

Inverse of D 

Speed of sound in the free stream 

r\j 

Constraining power of a constraint. If a = 1 
constraint is 100% effective. If a = 0 constraint 
has no effect at all 

Wing bending moment (about x-axis) 

Semi span of wing 

Aerodynamic coefficient (e.g. lift or moment coefficient). 

C is used as an experimental constraint 
e 

Lift coefficient 
Moment coefficient 

Wing root bending moment coefficient 
Reference chord length 
Section lift coefficient 
Section moment coefficient 
Section hinge moment coefficient 

Matrix that gives downwash (normalwash) in terms of pressures 
Deflection normal to lifting surface 
Unit vector in direction of axis 
Reduced frequency 


An average Mach Number between M 2 and M^ 
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wT 


x,y,z 


a 

P 

y 


Local Mach Number at wing surface 
Free stream Mach Number 


Average Mach Number between sending and receiving point 
Compressible radius (see Eqn. 73) 


Matrix that integrates pressures into aerodynami< 
parameters (e.g., C^> C^, c^, etc.) 

[S] [aC ] see Equation (9) 

P +• 


ci r a n r ».»n 


— «• /oa\ 


See Equation (45) 

[S] tgj-j-9 See Equation (54) 

See Equation (67) 

Weights given to the correction factors e for the 
minimization process, £ T e 2 = min 

Time 

Free stream velocity 


Local surface speed 

Correction factor = 1 + e (Called CF in computer program) 

Downwash (or normal wash) (Called W in computer program) 

Weights in the minimization process for estimates 

Cartesian coordinates right handed system 
x aft, y lateral (starboard), z vertical 

Angle-of-attack, also direction cosine for force or moment axis 


Direction cosine for force or moment axis 
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4C p 

aA 
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mod. 
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t 
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w 
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3/4 

1/4 


Dihedral of lifting surface 

Lower surface minus upper surface pressure coefficient 
Box area 

Incremental correction factors = W - 1 

Generalized incremental correction factors 

e = <{>e 

e /T 

Correction factor mode shapes 
Doublet potential function 
Circular Frequency 

Subscripts and Superscripts 

Stands for either p or w 

Designated or known correction factors 

Experimental 

Identifies estimates as opposed to constraints 
Hermetian transpose 
Modified values 

Identifies pressure modifying terms in the correction 
factor procedure 

Theoretical values 

Undesignated or unknown correction factors 

Identifies downwash modifying terms 

Deflection mode 1 

Deflection mode 2 

Three quarter chord point 

One quarter chord point 
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Matrix Notation 
Column Matrix 
Rectangular 
Diagonal 
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THEORETICAL DEVELOPMENT 


Basic Method 


Premultiplving Correction Factor Matrix . - A derivation of a real pre- 
multiplying correction matrix constrained to match aerodynamic data from a 
single downwash mode is presented in references 1 - 4. The use of Lagrange 
multipliers considerably simplifies the derivation so we present this alter- 
native derivation here. As an introduction we will rederive the same case 
first, i.e., the premultiplier for a single mode; then we will consider multi- 
ple modes and the postmultiplier. Whether the correction matrix is real or 
complex depends only on the experimental data: static data lead to a real 

matrix and oscillatory data lead to a complex matrix 

Assume that we have a matrix [A] of theoretical aerodynamic influence 
coefficients (AIC's) that relates the theoretical pressures {C p } on a set 
of aerodynamic finite elements to the dimensionless downwashes fw} at the 

k 

same aerodynamic elements by 


{ aC } = [A] {w} (1) 

p t 

The AIC's correspond to the reduced frequency of the experimental data and, 
hence, are real for static data and complex for oscillatory data. The pre- 
multiplying correction matrix [ W p J is used to obtain an estimate of the 

experimental pressure distribution {aC } from the theoretical distribution 

p e 

from 

uc > - r w p j uy (2) 

The subscript p refers to modification of the pressure distribution. The 
experimental force distribution is usually not known from the test data but 
only the integrated force and moment coefficients {C fi } are measured. An 
integration matrix [S] relates the experimental force distribution to the 
measured force coefficients through ; 
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{Cl = [s] (AC } 
e K e 


( 3 ) 


Combining equations (1) - (3) yields 


(C e > = [S] r W p J [A] (wJ 


( 4 ) 


which is the equation to be solved for the correction matrix T W p J given 

all the remaining terms in the equation. The remaining terms are all known: 

{C } and {w> are obtained from the test data, and [S] and [A] are known from 
6 

the mathematical model and the theoretical aerodynamic analysis of the con- 
figuration. In general, equation (4) is underdetermined, i.e., there are 
many more unknowns than equations. The method of least squares provides a 
solution. We require that changes in the theoretical load distribution shall 
be as uniform as possible or, in least-squares terminology, the weighted sum 
of the squares of the deviations shall be a minimum, where the deviation 
{e p > is defined as the difference between the correction factors -and unity. 

U p ) = (Hp-I) (5) 

We denote the weighting function by T p ; it will be discussed below. The 
weighted least-squares condition then becomes 


IT 


P 



V 


V V 


= a minimum 


( 6 ) 


where H denotes a Hermitian (complex conjugate) transpose. The Lagrange 
multipliers may be introduced by defining the error functional 

fp = 0/2) U p > H TT p j te p ) (7) 

and rewriting the measured generalized force coefficients (the constraints) as 
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{C e > = [S] ri ♦ Ep j {4C pt > 

= [S] UC p ) + [S] raCp^J {c p } 


The term [S] UC ) Is just the theoretical integrated pressures which are 
the theoretical cSefficients, tC t J . Thus 

{C f } = [S] UC ) < 8 > 

t p t 

Introducing also the following 

ty ■ [s] r v ( 9) 

. UC e ) = (C e J - tc t > 


gives finally 


(AC e > = [S p ] U p > 

The variation of the error functional f p is 

«fp * {c p )H fT p j Ue p ) 


( 10 ) 


( 11 ) 


and the variations of the incremental constraints, AC g , given by equation (10) 


are 


{«AC e > = [S p ] Ue p ) 
= 0 


( 12 ) 


The condition for the minimum subject to the constraints is then a linear 
combination of equations (11) and (12) set to zero in which the linear factors 
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are the Lagrange multipliers x p , 

6f + U1 H (6AC } = 0 (13) 

P P “ 

Substituting equations (11) and (12) into equation (13) yields 

<(«/ TT p J ♦ <»/ CS p ]) (Se p ) - 0 
Since the variation {6e p } is arbitrary, 

V" rT p J * {*/ [S p ] = 0 
or, after Hermitian transposition. 

fTpj <*p> + [S P ]H V = ° ( ' 4) 

The simultaneous soluti-on of equations (5), (10) and (14) yields the desired 
solution. The simultaneous solution leads first to the Lagrange k multipliers 
and then to e p as follows: 

<y - - (ts p ] rr p j-’ [S p ] H )-’ UC e } 

-fV' 1 [S P ]H ‘V 

and the correction factors are then 

(W p > = U>+ {e p > 

The premultiplying correction factors are written in a diagonal format 
use in subsequent aeroelastic analyses as in equation (2). 

The above results can be restated in summary form as follows: 

Solution of -*■ {AC g } = [S p ] {e p > 


(15) 

(16) 

(17) 
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2 

subject to ■* l e p Tp = min. 

(19) 

is 

- te) = [S p ] H ([S p ][S p ] H )'' (4C e ) 

(20) 

where 

[s p 3 - cs p ] 

(21) 

and 

( Ep j - 1 -rr p 

(22) 


The weighting function T p is arbitrary; the only requirement on it is 
that it should be positive. However, engineering judgment provides some 
guidance: if only a single constraint, e.g., the lift curve slope, is avail- 

able, one would prefer all the correction factors to be simply the ratio of 
its experimental value to its theoretical estimate. Accordingly, the recomm- 
ended choice for the weighting function for the premultiplier is 

{T p } = | [A] {I)| (23) 

However, other choices may be deserving of further investigation* 

Multiple Modes . - We next consider multiple experimental downwash modes. 
The derivation will be presented using two modes {w-j } and {W 2 }. For these 
two modes, the theoretical pressure distributions are 

{AC } = [A] {w-i } 

p tl 

{ aC } = [A] {w ? > 

p t2 

so the incremental experimental force coefficients aC^ and AC g ^ become 

uc el ) = <c e >- [S] uc Pt) j tvy 

= <C e >- [s p) ] <w p > ( 24 ) 
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a <c e >- [S p2 ] (W p ) 


(Continued) (24) 


These two equations may be combined into one set as follows: 

UC e ) = [Sp] (e p ) 

where 

UC e > 


and 

. ts P ] ‘ [s p,;V (26) 

k 

If again we impose the minimization condition 

T e 2 T = min (27) 

L P P 

then the solution is identical to equation (20) since equations (25) and 
(27) are identical to (18) and (19). 

Postmulti plying Correction Factor Matrix . - Now we consider the post- 
multiplying correction matrix. It is only necessary to consider a single 
downwash mode; the multiple mode case can be generalized by reference to 
equations (25) and (26). Since the postmultiplier modifies the downwash mode 
it defines an effective experimental downwash given by 

(Wg) = fW w J tw} 

in which the subscript w refers to modification of the downwash. Our new 
estimate of the experimental pressure distribution becomes : 


1 A S I 
= I I 



(28) 


{AC } = [A] {w } 

Pn e 


= [A] [W w ] (w) 

The experimental generalized forces or force coefficients are again given by 
equation (3) which with equation (28) becomes 



{C e J = [S] [A] n+e w J tw} 



= [S] {AC p ^} + ^ [A] fwj (e w > 

(29) 

And again noting 

that [S] {AC } = {C t > and introducing 

p t X 



1 — 1 
< 
l—l 

1 — I 
00 
1 1 

II 

I — 1 

00 

1 1 

(30) 

gives 

k 



{C e } - {C t } = [S w ] <£„) 


or since 

{C e J - {C t ) = {AC e > 



{AC e > - [S w ] (e w ) 

(31) 


Again the minimization condition is imposed, 

V e 2 T = min. (32) 

L w w 

The solution for e w is then identical to equation (20) since equations (31) 
and (32) are identical to equations (18) and (19). Again the correction 
factors are calculated from as follows: 

(W w > = {!> + (e w ) ( 33 ) 
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The weighting function T w is also arbitrary. However, the considerations 
that led to the recommendation of equation (23), in the premultiplying case, 
also lead to 


{T } * I U) T [A] | (30 

w 

which is to say that the weighting function is the lift coefficient induced 
by a unit downwash at each lifting element. Equation (34) is the recommended 
choice for the weighting function In the postmultiplying case, although other 
choices may still warrant further investigation. 

For multiple modes, say two, equations (31) and (32) provide the following 


<4 C e,> ” t S w, ] U w> 

{4C e 2 > = t S w 2 ] {£ w’ 

2 

V e T = min 
£ w w 

Again equations (35) and (36) can be combined into one as follows: 

(AC } = [S w ] (e w > 


(35) 

(36) 

(37) 


(38) 


where 


, 4C - 

l4C 4> ‘ . 1C . 


and where 


Equations (38) and (37) are now identical to equations (18) and (19) respect- 
ively and thus the solution is the same as before, i.e., equation (20). 
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Modifications to the Basic Method 


In some instances correction factors become unrealistic. In order to 
correct this situation when it occurs or to minimize the probability of its 
occurrence initially, various modifications can be introduced. Three such 
modifications are discussed here, i.e., estimates, correction factor modes 
and limits. Estimates are like constraints except that the "constraining 
power" can be varied. Correction factor modes constrain the distribution of 
correction factors such that the final distribution is a superposition of a 
limited set of well behaved, user input mode shapes. The "limit 1 feature 
constrains the correction factors (or any subset of them) to be above a given 
minimum and below a given maximum. 

Estimates . - In some instances data will be available in the form of 
estimates. These estimates can be based on past data, data from related 
configurations, two-dimensional data, empirical methods, or just past exper- 
ience. In any case they do not have equal weight with the experimental data 
considered so far. Consider the case where some experimental data are avail- 
able, leading to { AC g } , then the usual equation applies to these data: 

{AC e > = [S a ] {e a > 

where the subscript "a" stands for either p or w. If estimates exist, leading 
to {aC }, then it is desirable to minimize the difference between these esti- 
mates and the modified theoretical values. Let this difference be termed 
{e g }, then: 

(AC g ) « [S g ] {e a } + { £g > 

where {AC } = {C - C t > and S g is analogous to S fl with the exception that 

it refers^to the estimates C g and not the constraints C g . Thus we wish to 
minimize both (e a > and (e g > together and this is done as follows: 
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0 




This equation can then be solved in the usual manner producing the following 
result: 


I Te 2 + l = min 

If it is desired to give the e g values more or less weight in the minimization 
scheme, then the e g values must be weighted. 

I Te 2 + l (wT e g ) 2 = min. (40) 


The equations for the constraints then become 



(41) 


where 

U wT > - [wT] (Eg) (42) 


and where the values wT are the weights assigned to the errors e g . If the 
estimates are of high quality then the weights will be large. In the limit 
as wT ■* « , AC g becomes a constraint instead of an estimate and equation (41) 
reduces to the form of equation (18). Equation (41) can also be cast into 
the same form as equation (18) for the general case, i.e., wT finite, as 
follows: 


{ACg} = [S] (e) 


(43) 


with 
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{AC 


, . (i'll 

! IS I 


(44) 


[S] - 


s ! i_ 

9 1 wT 


(45) 


(e> 


■I 


1*A 

e WT ^ 


(46) 


Equation (40) can be written in terms of e as: 

i = min 


(47) 


where 


T for constraints 
1 


(wT)‘ 


for estimates 


Thus equation (43) and (47) are formally identical to equation (18) and (19) 
and thus have the same solution, i.e., equation (20). 


l ^ 

Currently the term is obtained from a term a where 


1_ 

wT 


1-a 

a 


(48) 


_4 

10 < a < 1.0 


where a is called the constraining power of the estimate C g . 
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Correction Factor Modes . - The correction factors can be expressed in 
terms of a set of modes <j> as follows: 

u> - M u> < 49 > 

Placing equation (49) into (43) gives 

uc e > = [S] [♦] {D 

= [S] { 1 } (50) 


where 


[S3 = [S3 [+] 


(51) 


If the minimization process is applied to e as usual 


I 


-2 

e 


min 


(52) 


Here the weight T is missing since it is usually not used with modes. Equa- 
tions (50) and (52) are then identical to equations (18) and (19) and thus 
have the same solution, i.e., equation (20). A similar expression exists for 
the postmulti plying correction factors. This approach allows a bias based 
on experience and past tests and physical reasoning, to be built into the 
correction factors. When estimates are considered equation (50) must be 
altered since the e wT are not fitted with correction factor modes. Thus 


{e} 



(53) 


and thus 


[S] 



(54) 
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and then the solution proceeds as before. 

Limits . - If certain basic properties of the weight factors are known, 
they could be limited to fall within a given set of bounds. If for instance 
the sign of an incremental weight factor is known to be positive, then it 
could be constrained to be positive. Also, for practical reasons, the maxi- 
mum value of the weight factors should be limited and thus the incremental 
weight factors are constrained to lie below this maximum. In general, the 
weight factors can be constrained to lie between a maximum and a minimum. 


e . < e < e 

min — - max 


(55) 


Notice that the generalized incremental correction factors, e, are the 
ones limited in the solution and not the actual ones, e. The values of e 
are the coefficients of the correction factor modes, <j>, and not the increment- 
al correction factors themselves. 

The basic procedure would be to set any generalized incremental weight 
factor to its maximum or minimum value if it exceeded these limits. This 
would require a multistep operation: (1) solving for the factors, (2) 

checking and setting those that exceeded the limits to the limit values, and 
(3) resolving. Before this can be done a capability must exist for assigning 
weight factors special values. This is easily accomplished as follows: 

«.> * <56) 

where [S] is defined in equation (51). The subscript u indicates those 
factors that are undesignated and d indicates those that are designated. This 
equation can be solved for {e u > in terms of the known quantities: 

<; d ) * [S B ] {e u > (57) 

Equation (57) effectively eliminates the designated factors from the minimi- 
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zation process. This equation can then be solved in the usual manner for 
{ } since {e^} is given. Specifically 


{ 



} = [S u ] l 'u } 


(58) 


where 

/ 4C \ ‘ <»<=,> - [S d ] (e d ) 
l mod > 


(59) 


The minimization scheme is then 


- * - 2 

T T e = min 
u u 


(60) 


Equations (58) and (60) are now formally identical to equations (18) and (19) 
and thus the solution is identical to equation (20). In the computer program 
the final e array that are modified or have reached their limits. is called 
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A New Postmultiplying Correction Factor Matrix 


The postmultiplying correction factor matrix developed in a previous 
section has been applied successfully to wings operating in pitch. Problems 
arise however when control surface modes are used. The discussion to be 
presented in the "Correlation Studies" section describes some of these prob- 
lems. As a result of these, a new postmultiplying correction factor matrix 
was developed and it is derived here. 

Viscous effects on airfoils can be thought of in terms of a displacement 
thickness added to the airfoil. The difference between the upper and lower 
surface displacement thicknesses produces a "decambering" of the airfoil or 
a change in the downwash w. 


The changes in downwash, fiw, exist over the entire airfoil or wing and not 
just in the region where w is non zero. These changes are a function of 
the pressure distribution on the airfoil. Usually the displacement thickness 
at a point is an integral function of the pressure distribution upstream of 
that point. If the general case of correction factor mode shapes is assumed 
then the downwash correction fiw can be expressed as: 

{6w} = U] {«e> ( 62 ) 

where {fie} is proportional to the integrated pressures , T^J. 

{fie} = { £ } (63) 

where [«J is given in terms of an integration matrix [N] and the pressures 
{AC p } . 

{*} = [N] {AC p } (64) 
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Combining equations (63) with (62) and using the result to obtain the 
corrected pressures leads to 


{AC } = [A] {w e > 

p e 

= [A] { w + [*] te» (65) 


The constraints {C } are obtained by integrating {AC } as follows: 
e 

(C p > = [S] {AC } 

e H e 

= [S] {AC } + [S *] {e } (66) 

p t v 


where 


[S p *] = [S] [A] C<t>] (67) 

Noting that {AC e ) = {C e J - {C t } equation (65) can be written as: 

(AC e > = [S p *] {e w > (68) 

Equation (68) has a form identical to that of equation (18) except [Sp] is 
replaced by [S *] and thus has the same solution, i.e., equation (20). Once 
found, {e w > can be placed into the expression for {AC p }, in equation (65), 
and the desired modified pressure found. 

Currently in the computer program the matrix [N] is simply either the 
identity matrix or the matrix [<t>] T . The identity matrix implies that the 
correction to the downwash is proportional to the local lifting pressure. In 
addition the above derivation is good for only one mode and thus the multiple 
mode option can not be used with the new postmultiplier. The new postmultiplier 
can be extended to multiple modes by simply replacing 

{ AC } withj__ e l( (69) 

AC e 2 
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and 


[S *] with 
P 


s * 

__r. 

S * 

L P 2 J 


( 70 ) 


but this has not yet been tried. 
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Transonic Effects Using Local Mach Number 


Empirical modification of theory is most meaningful if the theory quali- 
tatively matches experimental data. If the theory misses an important feature 
of the data the modified theory will also usually miss it. Transonic effects 
fall in this category. The classic lifting surface theory makes no provision 
for transonic effects and it is the purpose of this section to investigate 
some simple modifications to help remedy this situation. 

Direct Application of Local Mach Number . - Several methods based on the 
steady local Mach Number distribution have been tried and the results are 
discussed in later sections. 

One of these methods, discussed in the Introduction, consists of making 
a simple substitution of a local Mach Number distribution in place of its free 
stream value both in the kernel function and in the boundary conditions and 
pressure equation (see refs. 7 and 10). The local Mach Number distribution 
is taken from steady flow results. For the kernel function, application of 
a Mach Number distribution that lies somewhere in between the surface values 
and the free stream value was used. Tijdeman and Zwaan (ref. 10) suggest a 
local Mach Number distribution that lies half way in between the actual local 
and the free stream values. The reason for this is that acoustic signals 
propagate to the surface along various paths out in the fluid and thus propa- 
gate at some average between the surface value and free stream value. For 
the kernel, the local receiving point value of Mach Number and the free stream 
values were averaged and used in place of the free stream value. 

For the normalwash boundary condition and pressure evaluation the local 
Mach Number on the surface was used. Specifically, if M 2 U) is the local 
surface steady Mach Number distribution then the normalwash boundary condition 
w is (see ref. 10): 


w 

U 


M 2 (x) 




(71) 
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The second order Bernoulli Equation for steady flow is (see ref. 7) 


2 M ? (x) 

AC = AC [1 + p (-g 1)] (72) 

"so r CD 

where AC and AC are the first and second order pressures respectively. 

P Pso 

This method did not prove to be all that was hoped for and thus a second 
method was investigated. 


A New Transonic Effects Method . - In this section a derivation of the 
newly developed Douglas transonic effects method is presented. The basic 
method was conceived under the McDonnell -Douglas IRAD program however its 
implementation in two-dimensions and its application to the airfoil -control 
surface problem was done under the current contract. 


The lifting surface method is based on the following expression for the 
potential of a doublet, <j> d : . 


*d 


( e 1»t e 1X[H.(»-{)-R] I 
9n ( R I 


(73) 


where 


ojM 

00 

7u. 




2 2 
B r 


and n is the direction normal to the lifting surface. This expression can 
be rewritten as: 

>(t-x) } 

-J 


^d 3n ) 6 R 


(74) 


where 


T 


M 

oo 

B 2 U 


[R - M Jx-O] 


(75) 
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It can be shown that t has a physical significance. The term x is the 
acoustic time delay between the sending and receiving points. That is, x is 
the time it takes an acoustic signal, originating at the point S, n, c to 
reach the point x, y, z as the acoustic pulse washes downstream. 

This statement can be illustrated by the example of figure 1. 



Figure 1 


A signal is emitted at U, n, c) at time x = 0. At time x the wave front, 
traveling at the speed of sound, a M , has reached the receiving point at 
(x, y, z). During this time the wave center has travelled a distance U^x. 
Using the right triangle relations gives: 

r 2 ♦ U - x * U_ x) 2 - a 2 x 2 (76) 

Solving for x using the solution for a quadratic equation gives: 

T - \ ^ - X) + R) 

u/ 

which is exactly what is given in equation 75. 
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Thus t in the expression for <|> d has physical significance and it is the 
acoustic time delay between the sending and receiving points in a fluid moving 
with uniform velocity. This physical insight can form the basis of a correct- 
ion factor for the theory. For instance if the wave is in a flow field whose 
velocity varies in the longitudinal direction then the distance d travelled by 
the wave center is not U^t but is: 


■/' 


d = J U(t) dt ( 77 ) 

0 

If we consider U(t) to be made up of + 6U(t) then d also can be so split. 

d = U x + Sd ( 78 ) 


<5d = / 6U(t)dt 

J 0 


( 79 ) 


The wave center velocity U is being discussed, however this is not the velocity 
of the fluid particle located at the wave center as is the case for a 
uniform flow. The velocity U actually reflects the wave front speed and 
location and is the speed of an imagined wave center for the wave that strikes 
the receiving point. The wave front speed varies around the circumference of 
the wave, but the most important part of the wave is that part that strikes 
the receiving point. Thus the velocity U(t) is the time history of the wave 
center corresponding to that part of the wave that strikes the receiving 
point (x, y, z). As an approximation to the location of this part of the 
wave, it is assumed that it lies along a line connecting the receiving and 
sending points (shown dotted in figure 2). Thus 6U is the difference between 
the local velocity and the free stream velocity along the dotted line. If a 
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Figure 2 


coordinate R is defined lying along this line then the integral in time of 
equation (79) can be converted into a space integral in R as follpws: 

6d = f R 6U(R) 4dR 
J dR 

0 

dR 

where R is defined below equation (73) and where ^ is the speed with which 
the wave front moves along the radial coordinate. 

f = (..T 4 + U(R)D • i R 

The unit vectors i a , i and i R are defined in figure 2. 

For example, in the two-dimensional analysis for coplanar surfaces: 
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fid 



x 


= a - U(x) ( 8 °) 

r X M(x) - U^/a d - (81) 

J 1 - M(x) 

where M(x) is the local Mach Number distribution. As an approximation set 
U /a = U /a = M . 

oo' 00 00 00 

The time x can now be calculated using the right-triangle relations and 
the quadratic formula solution. 

[U^x + fid + K - x] 2 + r 2 = a 2 ? 2 

Solving for x gives: 

x = {M JK - x) + R} (82) 

U r 

oo 


dx 

dt 


Thus 


fid 


where 


(C - x) 2 + e 2 r 2 (83) 

= x - fid (84) 

It is immediately evident that x has exactly the same form as x (see 
equation (75)) except that x is replaced, in the expression for x, by x - fid 
in the expression for x. In essence then the receiving and sending points 
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have increased their separation (in the x-direction) as far as the acoustic 
time delay t is concerned. Is there any reason to carry this increase in 
distance to other parts of the potential function, specifically, to the radius 
term in the denominator? It seems so. First, this radial distance is 
already modified in the expression for x, see equation (82). Second it is 
known that transonic effects exist in steady flow (w = 0) where x has no 
effect; i.e., = |- (1/R). Thus it seems appropriate to add -fid to all 

x - £ terms. Thus 


J d (x - C, y - n, z - t, a), M^) = <t» d (x - 5 - fid, y - n* z - t, w, M w ) (85) 


where $ d is the potential modified for transonic effects. 

This method has been implemented for the two-dimensional case and the 
results are discussed subsequently. 


One variation of this method that is possible is to use an average Mach 
Number between sending and receiving points and define 6U as th& difference 
between the local value of velocity and this average. Thus the term M^ is 


replaced with M^ where 



1 



M(x)dx 


5 


(86) 


This method has also been tried and results using this variation are also 
discussed subsequently. 

A final consideration is the determination of the local Mach Number 
distribution, M(x). Tijdeman and Zwaan (ref. 10) note that the local surface 
Mach Number distribution should not be used but that some average between 
it and the free stream Mach Number should be used. This is because the 
signal arriving at a point has traveled both in the vicinity of the airfoil 
and out in the flow field. The recommendation of Tijdeman and Zwaan has 
been adopted in the present method. 
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Other work by Tijdeman and Bergh (ref. 11) can also be brought to bear 
on this work. Specifically a fully two-dimensional acoustic solution of a 
source pulse located at the control surface hinge was calculated for the 
case of a nonuniform flow field. This solution produced the exact time lag 
t from the hinge line to all other points on the airfoil. The equivalent 
distance x, and also R, from the hinge line to the receiving point can then 
be calculated using this information and the equation relating t to x. Thus 



If acoustic solutions were obtained for all other sending points then all 
the necesary x for this theory would be available. This method would be 
accurate, however it would require many expensive acoustic solutions. It 
appears that each of these solutions requires a computing effort comparable 
to a direct solution, by finite difference, of the original problem. This 
conclusion however, remains to be seen and further study is required. 
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CORRELATION STUDIES 


Local Mach Number Studies 

Several methods of accounting for local steady Mach Number variations 
in the oscillatory lifting surface theory have been studied for the two- 
dimensional case. A mathematical description of these methods has been 
presented in the Theoretical Development section. The cases considered 
here are for a two-dimensional symmetric airfoil (NACA 65A006) with an oscill- 
ating 25% chord flap. The local Mach Number variations over the airfoil at 
zero angle of attack are given in reference 11 . 

The first and simplest of the methods studied involves simply making a 
direct substitution of the steady local Mach Number in place of its free 
stream value. In general, this approach does not produce substantial changes 
in the pressures from their classical values. 

In figure 3 the symbols marked by triangles indicate the pressures 
calculated using the local Mach Number in the downwash boundary condition. 

The pressures are reduced from their classical values (indicated by dots and 
a dashed line) as expected, but not by very much. 

The circles indicate the pressures calculated using the local Mach 
Number in the downwash as well as in the kernel function. In the kernel 
function the average between the local receiving point Mach Number and the 
free stream value is used. The pressures again are generally reduced but not 
by any substantial amount. 

The use of local Mach Number in the second order Bernoulli equation 
(Ashley, ref. 7) also produces very little change. In figure 5 this change 
is observed as the difference between the circles and triangles. This change 
is about the same order of magnitude as the other changes except it is gener- 
ally in the opposite direction. 

The second method studied is new and is described in the Theoretical 
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Development section. The basic idea of this method is to transform the 
longitudinal distance between sending and receiving points depending on the 
time it takes an acoustic signal to travel that distance. A variation of 
this method is simply to replace the free stream Mach Number, M^, by M^, an 
integrated average of the Mach Number distribution between sending element 
and receiving point defined in equation (86). 

On the face of it the new method gives the best correlation when M m and 
not M is used. Figure 4 presents a comparison between the two methods at 
a Mach Number of 0.875. Near the leading edge the basic method designated 
"Present Method (M^)" and indicated by triangles produces the best agreement 
between experiment and calculation. Near the position of the steady shock 
wave however the peak pressure is better predicted by the variation of the 
basic method designated "Present Method (Mj and indicated by circles. The 
location of the calculated peak is slightly forward of the experimentally 
observed peak. Either method, however, is better than the classic theory 
(indicated by dots) for predicting pressures as comparisons with experimental 

data shown. 

Two features of the experimental pressure distribution illustrate trans- 
onic effects. One of these is the reduced leading edge pressure levels, and 
a second is the bump or peak in pressure near the location of the steady shock 
wave location. The new method qualitatively reproduces these features. 

However there is reason to suspect that the basic version of the new method 
(triangle) under-predicts the leading edge pressure. The reason for this 
lies in the fact that, even though the calculated and experimental pressures 
agree near the leading edge, viscous effects have not yet been accounted for 
and these effects reduce the calculated loads even further. A drop inthe 
leading edge loading caused by application of viscous effects to the (Mj 
variation of the basic method (circles), renders this method more acceptable 
than before. However these effects are not large enough to bring the calcu- 
lated pressures in line with the experimental values (see fig. 38). Further 
study is required in this area to decide which method is best or to discover 
other more accurate variations of the basic method. 
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The application of the new transonic method to a lower Mach Number, 

(0.85), is shown in figure 5. The agreement is good near the leading edge 
but only a slight indication of the shock bump is given by the theory. Also 
shown in this figure is the effect of local Mach Number on the Bernoulli 
equation (see equation (72)). The difference between the circles and triangles 
indicates this effect. 

All applications thus far have been for the steady case. Figure 6 
presents a comparison of the Present Method (new transonic theory) for the 
case of the control surface oscillating at a reduced frequency k^ = |jj = 0.059. 

oo 

Also shown in this figure is a calculation done using the Traci 6t al method 
(ref. 12) and a calculation done using the classic theory. The finite 
element theory of Traci et al predicts the bump at the shock wave fairly 
accurately however is not as good as the Present Method elsewhere. One part 
of the pressure distribution that does not seem to be predicted by any of 
the theories is the depth of the dip in pressure behind the shock. 

Figure 7 presents a comparison similar to that in Figure 6 *but at a 
lower Mach Number (0.85). Also instead of in ph ase (Real) and out of 
phase (Imaginary) parts given, amplitude AC p = y( ReAC p ) 2 + (ImAC p ) and 
phase angle = tan" 1 (ImAC p /ReAC p ) are presented. Again, as in the steady 
case the bump at the shock is barely noticeable in the Present Method and of 
course absent altogether in the classic theory. 

Figure 8 presents a comparison of the Present Method, classic theory 
and experimental data for a case similar to that presented in figure 7 
except that the Mach Number is 0.875 and the reduced frequency, k r = 0.176. 
Again pressure amplitude and phase angle are shown. The two variations of 
the Present Method are in better agreement with the experimentally obtained 
pressure amplitudes than is the classic theory. However the same can not be 
said of the phase angles. The Present Method follows the experimental phase 
angle curve from about the 40% chord on to the trailing edge. However none 
of the theories follows the curve forward of that point. Tijdeman and Bergh 
(ref. 11) present a modified phase angle curve based on a full two-dimensional 
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acoustic solution of a pulse located at the control surface leading edge 
(see section on New Transonic Effects Method). This approach gives very good 
agreement with the experimental phase angle data (see figure 31 of reference 
11). The correction was applied only to the phase angles and not the pressure 
amplitudes. The calculated phase angle was simply corrected using the 
additional time lag over and above that experienced in uniform flow. This 
additional phase lag was not used internal to the theory but applied after 
the theoretical calculation was completed. The section of this report 
entitled "A New Transonic Effects Method" describes how this acoustic type 
of information can be used internally with the theory so that the pressure 
amplitudes are also effected. This approach has not yet been tried. 

A possible explanation of the phase angle differences between theory 
and experiment might be due to the fact that signal fronts, which emanated 
from the control surface, do not exactly travel normal to the flow as assumed 
in the Present Method. Tijdeman and Bergh have shown that the wave fronts 
are actually inclinedto the flow to a considerable degree, within the super- 
sonic zone. This being the case the wave fronts impinge on the^ forward part 
of the airfoil (forward of the shock wave) with very little longitudinal time 
delay. This would explain the flattening of the phase angle curve in front 
of the shock. 

Thus far detailed pressure distributions have been discussed. Attention 
is now focused on the forces and moments these pressures produce. Figures 
9 through 12 present comparisons of the present method with the classic 
theory and experimental data. In figure 9 the Traci, Farr, Albano theory 
is also plotted. This figure shows that the Present Method (MJ is in better 
agreement with the data than is the classic method. As expected the Present 
Method and classic theory tend to coalesce at low Mach Numbers, out of the 
transonic region. The transonic peak lift, predicted by the Present Method, 
occurs earlier, as Mach Number is increased, than does the experimental data. 
Also the dip occurring after this peak is not nearly as deep as shown by 
the experimental data. 

Both theories show values of lift coefficient that are higher than the 
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experimental values. This is due to reduced flap effectivity caused by 
the viscous boundary layers. 

Figure 10 presents a comparison similar to the previous figure except 
that pitching and hinge moment coefficients are considered. Of particular 
note is the over prediction of hinge moment by both theories. Again this is 
due to viscous effects on the flap. 

The last two figures have dealt with force and moment coefficients in 
steady flow for an airfoil with a deflected flap. Figures 11 and 12 present 
the same data for the oscillatory case. (The reduced frequency varies from 
0 098 at M = 0.5 to 0.057 at M = 0.901 in these figures.) Generally speaking 

the same trends and conclusions hold for these figures as for the previous 
two figures. 
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Subsonic Cases 


Oscillating Ming with Control Surface . - Extensive low speed wind 
tunnel measurements of static and oscillatory pressure data have been made by 
Hertrich (refs. 14 and 15) on straight and swept wings with a full span con- 
trol surface. The wings had no taper and the control surface had a 30% chord; 
two aspect ratios, 2.5 and 3.1, were tested by changing the exposed span in 
the tunnel. The swept wings had a sweep angle of 25 . A later oscillatory 
test of the swept wing was made by Forshing, Triebstein, and Wagener (ref. 

16) in which the full span control surface was split approximately in half 
(the inboard flap had 46.59% of the span( and the aspect ratio was set at 

2.94. 

The pressure data from the first tests (refs. 14 and 15) were integrated 
by Hertrich to obtain lift and moment coefficients and the static values for 
the swept wing with aspect ratio 3.1 have been used here as constraints to 
determine correction factors. The static values corrected for wind-tunnel 
wall interference are:' lift curve slope C L = L/qs = 3.13 per radian, 

pitching moment curve slope C m = M/qSca = 0.148 per radian, flap lift effect- 
iveness C, = L/qS6 = 1.95 per a radian, flap pitching moment effectiveness 

C = M/qSc6 = -0.432 per radian, and the flap hinge moment coefficient 

m. 

C h = H/qSc6 = -0.0350 cos 25° = -0.03172 per radian where the sweep correct- 
ion is added to give the moment about the hinge axis; the reference area is 
S = 0.564 m^, the reference chord is c = 0.6 m, and the pitch axis is located 
at 61.5% of the root chord. 

Rolling moment coefficients were not derived from the pressure data and 

neither was the hinge moment due to angle of attack C h . The available data 

a 

permitted a maximum of five constraints, and two sets of constraints were 
investigated; the first set used two constraints from the angle of attack 
data, C, and C , and the second set used all five constraints. The use 

of the flap rotation data alone without the angle of attack data was not 
considered. 
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The theoretical basis for the correction factors is the Doublet-Lattice 
Method (DLM) of reference 17. The idealization of the lifting surface consists 
of 110 boxes resulting from 11 strips and 10 equal chordwise divisions. The 
chordwise division on each strip result in 7 boxes on the primary surface and 
3 boxes on the control surface. The strip widths ay i are chosen so that the 
strip centerlines fall along the lines of pressure taps. The span of 0.940 m 
is divided into the following strip widths from root to tip: Ay ] = 0.110 m, 

ay 2 = 0.080 m, ay 3 = 0.075 m, ay 4 = ay 5 = ay 6 = Ay ? = ay g = ay 10 = 0.090 m, 
and Ay n = 0.045 m. The pressure stations correspond to the strips as 
follows: pressure station VII is on Strip 1 , VI on Strip 4, V on 6, IV on 8, 

III on 9, and finally station II is on Strip 10. Pressure station I is too 
close to the tip to permit a meaningful calculation. 

The theoretical pressure distributions are compared to the experimental 

measurements In figures 13 through 24. The theoretical estimates of the five 

constraint parameters are: C. * 3.207462. C = 0.179494, C. ■ 2.131577, 

• a a 6 , 

C = -0.463554, and C h = -0.057784. Three additional parameters are also 

of 6 interest. These are 6 the locations of the spanwise aerodynamic centers for 

angle of attack, y fl /s, and for flap deflection, y fi /s, and the hinge moment 

coefficient for angle of attack, C h . Their theoretical estimates are 

a 

y /s = 0.452071, y./s = 0.464614 and C. = -0.021034. 

''a o a 

A typical set of correction factors is shown in table I; it is for a 
premultiplying matrix and is based on five constraints. The factors are 
listed in order from leading edge to trailing edge on each strip beginning at 
the root; factors 1 to 10 are on Strip 1, and factors 101 to 110 are on Strip 
11. The first seven factors on each strip apply to the primary surface and 
the last three apply to the control surface. The general trends seen in 
table I are a spanwise increase in factors from root to tip and a chordwise 
increase toward the hinge line. The minimum correction factor in table I 
is a 0.255873 for Box No. 10 and the maximum factor is 2.00893 for Box No. 108. 
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TABLE I 


^rppcrjlCN r aCTO^S * * PPE^LLTIPLIEP CASE 


1 

0 . 84 8 2 78E + 00 

C.O 

2 

0. 6639 80E *00 

0.0 

3 

0.877532E+00 

0.0 

A 

0.88593 5E+CC 

C.O 

5 

C.884442E+CC 

0.0 

6 

0 . 863 395E + CC 

0.0 

7 

0.795167E+0C 

C.C 

e 

0.35126 IE + C0 

0.0 

9 

0.377085E+C0 

0.0 

LO 

0 . 2 55 8 73E ♦ CO 

c.c 

1 1 

0 . 8 1 706 5E+0C 

0.0 

12 

0 • 894 144E + CO 

0.0 

13 

0.907865E+0C 

C.C 

1 4 

0.516355E+GO 

0.0 

15 

0.916132E+03 

0.0 

16 

C. 899734E+00 

C.O 

17 

C.845436E+00 

0.0 

18 

0.431832E+00 

0.0 

19 

0.430533E+CC 

c.c 

20 

C . 2 5 5 866E + CC 

0.0 

21 

0.90 l 404E ♦ 00 

0.3 

22 

C.920260E+0C 

c.c 

23 

0.935662E+00 

0.0 

24 

0.946 28 7E +03 

0.3 

25 

0. 94949 8E + C C 

c.c 

26 

C .935445E+C0 

0.0 

27 

0.897291E+00 

0 .0 

28 

0.510927E+00 

C.O 

25 

C.475718E+00 

0.0 

30 

0. 330467E+00 

0.0 

31 

0.928828E+OC 

C.O 

3 2 

C.95CCG3E+C0 

0.0 

33 

0.967955E+ 30 

0.3 

3A 

0.981 8C3E+C0 

c.c 

35 

C.S856C2E+C0 

0.0 

36 

0. 98678 7E+00 

0.0 

37 

0.95819 3 C +C0 

c.c 

38 

C.6C7972E+C0 

0.0 

39 

0.530440E+00 

0.0 

AO 

O.366155E+0C 

c.c 

A 1 

C . 5 6C 566E + CC 

0.0 

42 

0 .984268E+03 

0.0 

A3 

0. 13052JE+J1 

C.O 

4 A 

C . 1C2273E+C1 

0.0 

45 

0. 103556E+01 

0.0 

A6 

0. 104067E+01 

C.O 

A7 

C.102777E+CI 

0.0 

48 

0.721815E+00 

0.0 

A9 

0.591 156E + CC 

C.O 

5C 

C.4CC92CE+00 

0.0 

51 

C.993722E+33 

0.0 

52 

J. 1 ) 15 83E + C l 

C.O 

53 

C.104356E+C1 

0.0 

54 

0. 106452E+C1 

0.0 

55 

0. 1082C9E+01 

c.c 

56 

0.1C9505E+01 

0.0 

57 

0. 109904E+01 

0.0 

58 

C.840227E+CC 

C.O 

59 

C.6467C1E+00 

O.C 

60 

0.424058E+0C 

0.0 

61 

0. 102303E+01 

C.O 

62 

C . 105634E+JI 

0.0 

63 

0. 1G8254E+01 

0.0 

6A 

0.1 10650E+C 1 

C .0 

65 

0. 11284 7E +01 

0.0 

66 

0.11495RE+01 

0.0 

67 

0.1 17262E+01 

C.O 

68 

0.96 558CE + CC 

0.0 

69 

0.692697P+03 

3 . 0 

7J 

0.428853E+00 

0.0 

71 

C . 106328E+01 

0.0 

72 

0. 109351E+CI 

0.0 

73 

0.1121 73E + C 1 

O.C- 

7 A 

0.1148 14E + 01 

0. 0 

75 

0.117424E+CI 

0.0 

76 

0. 120438E+C1 

C.O 

77 

C.125149F+CI 

0.0 

78 

C. 1 10761E+01 

0.0 

79 

0.726251E+0C 

0.0 

8C 

C .4 1C 725E + CC 

0.0 

81 

0. 109930E+01 

0.0 

8? 

0.1 13103E+01 

c.c 

83 

0.1 16C48E+0 l 

0.0 

84 

0.11 8864E+0 l 

0.0 

85 

0. 121911E+01 

O.C 

86 

C. 126123E + C1 

C.O 

87 

0.134366E+01 

0 .0 

88 

0. 129296E+C l 

0. J 

85 

C.745826E+00 
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0.126399E+01 

0. 0 

96 

0.132572E+C1 

0.0 

97 

0. 146760E+31 
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C.O 
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0.0 
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0.119185E+01 

0.0 
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0.121559E+C1 

C.O 
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C. 125326E+01 

0.0 

105 

0.130061E+CI 

0.0 

106 

109 

0. 1386C9E+01 
0. 945436E + C0 

C.O 

0.0 

107 

110 

0. 159604E+01 
0.6153065+00 

0.0 

0.0 

108 

3.230893F+C1 

0.0 




TABLE II 



Type of 

Number 

V s 

V s C h 

Ch 

h 6 

Correction 

of 


ct 


Matrix 
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None 

0 

0.452071 

0.464614 -0.021034 

-0.057784 

Pre- 

2 

0.456751 

0.469814 -0.021746 

-0.059959 

Pre- 

5 

0.484939 

0.522318 -0.010117 

-0.031721 

Post- 

2 

0.453096 

0.465632 -0.022409 

-0.057711 

Post- 

5 

0.469176 

0.491720 +0.013014 

-0.031721 
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The modified static pressure distributions for angle of attack are shown 
in figures 13 and 14, and for flap deflection are in figures 15 and 16. 

Perusal of these figures indicates the following results. For the angle of 
attack loading, both the premultiplying and postmultiplying corrections move 
the theoretical results slightly away from the experimental data, the post- 
multiplier causing a little greater change. The effect of five constraints 
is greater than that of two. For the flap loading, both the pre- and post- 
multipliers based on two constraints have small effect. The corrections 
based on five constraints improve the correlations on the control surface but 
increase the discrepancies on the wing. The postmultiplier causes a much 
larger change and, although the data show a pressure reversal near the trail- 
ing edge, the postmultiplier exaggerates this reversal to the extent that the 
sign of the hinge moment is reversed. Table II shows the effects of the 
four correction matrices on the aerodynamic centers and hinge moments. All 
of the correction matrices resulted in an outboard shift of the aerodynamic 
centers, the largest shift coming from the 5-constraint premultiplier. 

Two constraints did not improve the hinge moment prediction, and the 5- 
constraint postmultiplier lead to an unreasonable prediction of C ha . The 
effect of additional constraints based on estimates is a topic deserving 
further investigation. 

We can anticipate similar discrepancies when the correction factors 
derived from static data are applied to the oscillatory cases, and, indeed, 
they are shown in figures 17 through 20 for the angle of attack oscillating 
at k = 0.622, and in figures 21 through 24 for the flap oscillating at 
k =*" 0 . 752 . The theoretical loading for the oscillating angle of attack is 
not changed significantly by either the premultiplier or the postmultiplier 
based on two constraints and both the real and imaginary parts are affected 
about the same. In some regions the theory is shifted toward the data and in 
others the theory is moved away from the data. The effects of five constraints 
are more extreme. The 5-constraint premultiplier improves the correlation 
for the real part but only improves the agreement for the imaginary part on 
the control surface while diverging on the wing. The 5-constraint postmultiplier 
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is substantially worse in correcting the real part but is no worse than the 
premultiplier in modifying the imaginary part. Again, the theoretical load 
distribution from the oscillating control surface is not changed significantly 
by either the 2-constraint pre- or postmultipliers. However, some improvement 
is noted with the 5-constraint corrections although it is only slight. As 
in the static case, an outboard shift in loading occurs with all correction 
matrices and for both modes of motion. 

The above applications of correction matrices have achieved very limited, 
if any, success. The lack of improvement in the most elementary case, however, 
is rather puzzling. This was the case of the static loading at angle of attack 
for which the correction factors were derived using the two constraints of 
lift and pitching moment. The pitching moment constraint was expected to 
shift the theoretical chordwise center of pressure in such a manner that the 
predicted pressure distribution would be closer to the experimental data. Two 
explanations for the lack of improvement appear possible. The first is that 
the theoretical loading. in the leading edge region differs so much from the 
data that it dominates the correction factor calculation and results in a 
distorted loading. The second possibility is that the limited number of 
pressure taps near the leading edge prevented an accurate evaluation of the 
leading edge contribution to the pitching moment. A strain gage measurement 
of pitching moment would have shed some light on this possible difference. 

A number of options were not pursued with these data which may have shoy/n 
better correlation. First, only one configuration was studied here, the 
swept wing with aspect ratio 3.1; as noted above, straight wings with two 
aspect ratios and a swept wing with another aspect ratio were also tested. 

Next, only the reported integrated loads were used as constraints: the two 

angle of attack coefficients and the three additional control surface coeff- 
icients. The three control surface coefficients were not used as constraints 
by themselves, nor were additional constraints used based on estimates of 
rolling or bending moments. The new postmulti plying matrix was also not 
investigated. Finally, it would have been interesting to apply complex 
correction factors derived from the oscillatory angle of attack data at 
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k p = 0.622 to the oscillatory control surface data at k r = 0.752; however, 
this would have required integration of the published oscillatory pressure 
data to obtain the complex constraint coefficients. 


Ming With and Without Leading Edge Droop . - Trailing edge control surfaces 
are studied in several other sections of this report. In this section, an 
attempt is made to study leading edge control surfaces. Usable data for such 
devices is very scarce. Several references have been investigated; however, 
only reference 18 proved in any way useful. The leading edge device described 
in this reference is. a wing droop of 6°. The droop was applied to the first 
19% of the wing chord along its entire span (see fig. 25). The idealization 
shown in this figure is for the Doublet Lattice Method (DLM). The fuselage 
was simplified as simply a wing extension to the centerline. 


A steady case at M = 0.80 is considered and the uncorrected calculated 
results using the DLM for a = 4°, (no droop) agree very well with the experi- 
mental data (see fig.. 25). Only a lateral shift in the center of pressure 
seems evident. Correction factors were developed for this case to correct 
this slight deviation in the theory. The constraints used are lift, pitching 
moment and bending moment coefficients. These coefficients were summed on 
strips outboard of the station y/(b/2) =0.11 and are defined as: 
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A premultiplier and a postmultiplier (new type) were tried with equally 
good results on the span loading. Figure 26 illustrates the effect of the 
correction factors on the spanwise distribution of aerodynamic center. The 
correction factors increased the accuracy of the aerodynamic center inboard, 
but decreased it outboard. 

The single mode application of both pre- and postmultipliers, also 
produces good results for leading edge droop span loadings (fig. 27). Notice 
that the unmodified results are approximately half of the experimental values. 
The experimental data were difficult to read on the plots (open squares). 

Thus, the pressure distributions were integrated to produce the darkened 
squares. The pressure distributions themselves were difficult to” integrate 
accurately since there were down loads at the nose and uploads near the bend 
in the chord, such that the total loads were small. If the correction factors 
possessed only a slight variation in the chordwise direction, the balance of 
integrated load could shift drastically as a percent of the total. 

The flow field near the wing changes at approximately a = 8°. Here the 
flow is to a large extent separated from the upper outboard surface. A 
comparison of uncorrected theory and experimental data, for the case of 
a = 10°, no wing droop, in figure 28, shows a loss of lift outboard of the 
40% semi -span. Application of both pre- and postmultipliers (New), using 
C^, and Cg (bending moment at the centerline), show a much improved pre- 
diction of span loading. Also, shown in this figure is the application of the 
premultiplier correction factors, obtained at a = 4° to the a = 10° case 
(diamond symbols). The span loads are improved which shows that data obtained 
at one angle of attack can be profitably applied to other angles of attack. 

The corrections generated at a = 4^ are not as large as those generated at 
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a = 10° because flow separation exist in the latter case. However, both . 
corrections are in the same direction. Therefore, application of correction 
factors for a = 4° improves the results for the 10° case. In general, the 
reverse may not hold; i.e., the correction factors obtained at a = 10 (or 
larger angles) may be too large and an excessive correction may result leaving 
the corrected data further from the experimental data than there were originally. 
It does seem safe, however, to apply correction factors obtained at one angle 
of attack to other nearby angles if the flow is qualitatively similar (e.g., 
no great changes in flow pattern). 
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Transonic Cases 


In this section applications of the correction factor technique are 
made to the same cases considered in the "Local Mach Number Studies" section. 
Specifically a two-dimensional symmetric airfoil (NACA 65A006) with an oscill- 
ating 25% chord flap is used. 


Figure 29 illustrates the difference in results obtained when the classic 
theory (subsonic compressible) and the new transonic theory (Present Method 
(M )) are corrected. A premultiplying set of correction factors were obtained 
using three constraints; lift, moment (1/4 chord) and control surface hinge- 
moment (3/4 chord). Each theory was corrected to the proper experimental 
constraints, i.e.. 


m 


1/4 


4.93 

-1.57 


M = .875 

oo 

k r = 0.0 


c. = -0.053 

h 3/4 

where the characteristic length is the chord and the downwash over the control 
surface is unity. The classic theory does not have the bulge in pressure, 
near the compression (or shock) region for the steady flow as does the experi- 
mental data and applying correction factors will not make it appear. Thus 
correction factors can not make a qualitative feature appear where none existed 
before. The corrected classic theory does not compare well with the experi- 
mental data and the correction factors themselves, (1 + e), show fairly 
large deviations from unity especially near the leading and trailing edges. 

The Present Method (M^) however possesses a qualitative similarity with 
the experimental data and thus it is a better candidate for correction. 

Figure 29 shows such a correction. The bulge in pressure as calculated by 
the present method is amplified as it should be. The loading on the flap 
however is reduced, again as it should be, however the shape of the flap load 
is distorted. The correction factors themselves are better behaved for the 
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Present Method, (M ), showing large deviations from unity only on the flap 
surface. 

Figure 30 presents the results of applying the premultiplying correction 
factors, obtained for the steady case, to an unsteady case. That is, the 
correction factors shown at the bottom of figure 29 are applied to the 
oscillatory results of the Present Method (M^) and the classic theory 
(k = ojc/ 21^ = 0.059). Since the correction factors are real they do not 
effect the phase angles of the pressures but only the amplitudes, |AC p |. Also 
shown in the figure is a pressure distribution corrected using factors based 
on the complex lift moment (1/4 chord) and hinge moment obtained for the 
unsteady case. Specifically 


C £ 

3.5 - i 1.18 

M 

00 

= 0.875 

Cm l/4 

-1.66 + i 0.07 

k r 
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-0.057 - i 0.016 



h 3/4 





The correction factors obtained in this manner produce pressures ‘that are 
close to those obtained using the steady correction factors (except near the 
flap) even though the constraints in lift are considerably different in the 
two cases. There is one slight anomaly in the phase angle for the complex 
constraint case (k f = 0.059) and it exists on the last two pressure points on 
the flap. The phase angle there is quite large however these angles do not 
have a large effect since the amplitude of pressures is very small there. 

The question arises; to what extent can static correction factors be 
applied with accuracy to the oscillatory case? Figures 31 and 32 illustrate 
the effect of static correction factors on lift, moment and hinge moment 
coefficient versus reduced frequency for a Mach Number of 0.85. Considering 
first the lift coefficient it is noticed that the accuracy of the imaginary 
part is increased up to k f ■ 0.2. Beyond this point application of correction 
factors decreases the accuracy of the theory. For the real part of the lift 
the corrected theory is more accurate onlybeldw a reduced frequency of 0.06. 
For the pitching moment and hinge moment the cross over point is roughly 
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k =0.1. On the average then the static correction factors are useful up. 
to about a reduced frequency of 0.1. Beyond this point it is better to use 
the original theory. 

It is probably true that the accuracy of extrapolating correction factors 
versus reduced frequency depends on Mach Number. Figure 33 gives an indication 
that as Mach Number is reduced the accuracy increases. Specifically, static 
correction factors have been applied to the oscillatory case (k r = 0.098) with 
very good results. Both amplitude and phase angle are improved. 

The theory used in figures 31, 32 and 33 is a variation of the new 
transonic method presented previously. Specifically the variation utilizes 
an average local Mach Number (MJ in place of the free stream Mach Number, 

M . Fiqures 31 and 32 show the application of two separate types of correct- 
ion factors; a premultiplier type, the type used in figures 29 and 30, and a 
postmultiplier type. The postmultiplier is actually an additive viscous type 
of correction. It can be seen that this correction does not extrapolate to 
higher frequencies as well as the premultiplier type (as far as the lift 
coefficient is concerned). As the frequency is increased the experimental 
data approach the unmodified theory. One interpretation of this fact is that 
as the frequency is increased viscous effects are reduced. 

The postmulti plying correction factor (designated as "New Post") used 
in figures 31 and 32 is the new postmultiplying correction factor discussed 
in the Theoretical Development section. The reason a new type of postmultiplier 
was needed is because the original one seemed to fail when control surfaces are 
considered. Postmultipliers correct the downwash matrix. When all downwash 
values are non zero, e.g., wing pitch, the method seems to work. However, 
when this is not the case, e.g., control surface deflections, the method 
fails entirely. The corrected downwash values are either large and erratic 
themselves or they cause large and erratic pressures due to the modified 
downwash. 

It was hoped that the introduction of correction factor mode shapes would 
smooth out the corrected downwash and produce accurate results. This did not 
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happen. Even though smooth, well behaved functions were used the results were 
unrealistic. Although not tried, it seems that limiting the maximum and minimum 
values of the correction factors probably would not help very much either. 


This failure of the postmultiplying correction factors led to an interest- 
ing investigation and subsequent development of the "New Postmultiplier . 

The investigation consisted in finding out what downwash in the theory would 
produce the experimental pressure distribution. Specifically the theoretical 
influence coefficient matrix was multiplied by the vector of experimental 
pressures to produce a vector of downwash values. 


Figure 34 shows the results of this type of analysis (designated as 
Experimental) for a steady subsonic case (M^ = 0.5). Also shown is the 
theoretical downwash, i.e., unity over the flap. One thing is noticed immed- 
iately, there is a change in downwash ahead of the flap even though it is 
theoretically zero there. This downwash change is like a negative pitch of 
the entire airfoil. Figure 35 shows the camber (designated = 0.5) associat- 
ed with the downwash given in figure 34 and indeed it is like a negative or 
nose down pitch. This fact suggests that an additive type of correction 
factor, whereby all downwash values are changed, is necessary. This resulted 
in the development of the "New Postmultiplier" as described in the Theoretical 
Development section. This name is somewhat of a misnomer since the correction 
factor is additive and not multiplicative although the correction factors 
are proportional to the theoretical pressures. 


The results of applying the new postmultiplying correction factors are 
also shown in figure 34. Again lift, moment (c/4) and hinge moment (c3/4) 
coefficient were used as the constraints 


"m 


'l/4 

*3/4 


= -.70 


= -.0528 


The corrected values of downwash (circular symbols) agree well with the experi- 
mentally deduced downwash. The disagreements at the leading edge of the air- 


51 


foil and ahead of the flap are due to the fact that downwash is a sensitive 
function of pressure and slight variations cause large variations in downwash. 
With this in mind the agreement is very good especially over the flap itself. 

Applying this corrected downwash to the theory produces the results given 
in figure 36 for the pressure distribution. The results of the New Postmulti- 
plier agree very well with the experimental pressures. For reference, correct- 
ions by a premultiplier are also shown and these are also very good. The 
uncorrected theory is also presented for reference. 

At the low Mach Numbers used in the last few figures (M^ = 0.5) transonic 
effects are not present and any differences between theory and experiment are, 
in all probability, due to viscous effects. Figure 35 has shown that viscous 
effects modify not only the downwash over the flap but also over the forward 
part of the airfoil as well. This comes about due to the fact that the de- 
flected flap causes an induced upwash over the forward portion of the airfoil 

which in turn generates, a difference in boundary layer displacement thickness 
on the upper and lower surfaces. This difference in displacement 1 thicknesses 
causes an effective nose down pitch. 

It stands to reason that the correction factors generated at M^ = 0.5 
could be used to increase the accuracy of the theory at all Mach Numbers since 
viscous effects are present at all Mach Numbers. Figure 35 shows the effective 

cambers at M^ = 0.5 and M^ = 0.875 using the (Mj variation of the new transonic 

theory. Notice that the transonic camber can be thought of as composed of two 
pieces; one viscous piece very similar to that found at M^ = 0.5 and one trans- 
onic piece with the shape of a bump. This indicates that the accuracy of the 
corrected camber (or downwash) at transonic Mach Numbers can be increased if 
the subsonic (M^ = 0.5) results are known and used since it represents one 
part of the correction. 

Figures 37, 38, 39 and 40 give examples of applying correction factors 
obtained at M^ = 0.5 to other Mach Numbers for both pressures and aerodynamic 
coefficients. Specifically figures 37 and 38 present the results for Mach 
Numbers of 0.85 and 0.875 respectively. Up to three separate corrected 
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pressure distributions are shown in each figure. One is the result of applying 
a premultiplying correction factor matrix to the (Mj variation of the new 
transonic method. A second is the result of applying a new postmultiplier to 
the same theory; and third is the result of applying a new postmultiplier to 
the (Mj variation of the new transonic method. The last pressure distribution 
is seen to be the most accurate and a definite improvement over the unmodified 
(MJ theory (see fig. 4). 

Figures 39 and 40 give a clear picture of the effect of applying 
corrections obtained at M = 0.5 to other Mach Numbers. Figure 39 presents 

oo 

the lift coefficients associated with corrected and uncorrected pressure distri- 
butions. Two types of corrections are used; both pre- and postmultiplier 
(New). The theory used is the (Mj variation of the new transonic method. 

Figure 40 presents similar results for the pitching moment and hinge moment 
coefficients. The corrections developed at M^ = 0.5 greatly improve the theory 
as far as the lift coefficent is concerned. The pitching moment is not changed 
much because it was very close to the data to begin with. The hinge moment 
also is not changed much. 

Figures 39 and 40 show that corrections obtained at low Mach Numbers 
can be applied to the theories to improve accuracy at higher Mach Numbers. 

Figure 41 presents the results of correcting the theory with both a post- 
multiplier (New) and a premultiplier. First the theory is corrected using a 
new postmultiplier obtained at = 0.5. This represents a viscous type of 
correction. A premultiplier is then applied to the previously corrected 
results to account for transonic effects. This process produces a pressure 
distribution that approaches the data more closely than any of the others when 
it is combined with the (MJ variation of the new transonic method. 

Figure 42 presents typical correction factors for the steady two-dimen- 
sional cases considered in this section. The theory used is the Present 
Method (M ). There is a greater change in premultiplying correction factors 
between M = 0.85 and 0.875 than there is between M = 0.5 and 0.85. 

oo 
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Supersonic Case 


The arrow wing, shown in figure 43, has been chosen to illustrate the 
application of the correction factor technique to the supersonic case. The 
Douglas Supersonic Doublet Method (SDM) (ref. 20) has been used to determine the 
theoretical loads. The box idealization used is also shown in the figure. 

Notice that the tip of the wing has been clipped to reduce the number of boxes. 
Two modes are considered; (1) pitch (a = 4°) and (2) camber. The wing is 
operating at a Mach Number of 2.05 and a reduced frequency of zero. 

Figure 43 presents a comparison of uncorrected theory (dotted line), 
corrected theory and experimental data. The experimental values lie below 
the theoretical (uncorrected) values over the entire span. 


Four different methods of correcting the theory were tried. A pre- and 
postmultiplier (New) were applied using the pitch mode only and the results are 
very encouraging. The only real difference between the theory, corrected in 
this manner, and the experimental data appears at the wing tip. multiple 
mode case was tried using the pitch and camber modes and the results are good 
but not as good as the previous two corrections. The fourth method is the 
application, to the pitch case, of a premultiplier correction factor matrix 
that was derived for the camber case. Thus a correction factor derived for 
one mode (camber) is applied to another mode (pitch). The results are not 
very accurate on the inboard part of the wing but agree as well as the other 
methods on the outboard part. The constraints used are summarized as follows: 
Pitch a = 4° Camber 
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A similar set of corrections were applied to the camber case and the 
results are shown in figure 44. Specifically pre- and postmul tiplying correct- 
ion factors were obtained using the camber mode. In addition a multiple 
mode (pitch, camber) premultiplying correction factor matrix was derived 
using six constraints; i.e., C L , C M , C g for both modes. All three of these 
corrections give approximately the same good results except right at the 
wing tip. 

Two other types of correction factors are applied to the theory and these 
refer to applying the correction factors derived for pitch to the camber mode. 

On the inboard portion of the wing these correction factors over-correct the 
theory, but are accurate on the outboard portion of the wing. On the inboard 
portions of the wing the correction factors move the corrected theory further 
from the data than itwas originally in its uncorrected state. 

Figure 45 illustrates how correction factors modify the pressure distri- 
bution over the wing (in pitch) using a premultiplier. It seems that the 
reduction in lift due to the correction factors is taken out at the trailing 
edge rather than the leading edge as the experimental pressures would indicate. 
This fact might be explained if the experimental pitching moment were inaccurate. 

The postmultiplier does not directly modify the pressures but modifies 
the downwash. Modified downwash can be expressed in terms of modified camber. 

It is of interest to know how the postmultiplier (New) modifies the wing 
camber and figure 46 presents such a modification. The camber is reduced, by 
the correction factors, over most of the wing just as expected, since the 
boundary layer and separation regions act to reduce the effective wing camber. 

The postmultiplier, then, acts in a way that is consistent with physical 
processes. 
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RECOMMENDATIONS FOR DATA ACQUISITION 

Most of the data utilized in this study were pressure data, and the 
correction factors were derived using constraints that were obtained in some 
instances from integrations of the pressures. Certain errors are associated 
with integrations of pressures to obtain generalized forces, arising primarily 
from the limited number of pressure pickup points on a practical model. The 
forces and moments should be measured directly in addition to the pressures. 
Control surface hinge moments should be measured and so should rolling moments, 
i.e., root bending moments, because of the importance of the spanwise aero- 
dynamic center location. For swept surfaces it would be desirable to measure 
not only pitching and rolling moments in the streamwise coordinate system but 
also root bending moment and torque about some swept coordinate system, e.g., 
the 25% or 50% chord lines. 

Two significant deficiencies were observed in available experimental 
data besides the absence of combined pressure and force data. One was a lack 
of any systematic variation in reduced frequency in covering the range from 
steady flow to high frequency, i.e., k p of order unity. The correction factors 
are frequency dependent, and it is not reliable to use factors derived from 
low frequency data to predict pressure distributions at high frequencies. The 
second deficient area is the effect of Reynolds number. An important source 
of discrepancy between theory and test is the neglect of viscosity in the 
theory. When extensions of oscillatory lifting surface theory are made to 
account for viscous effects, data will be needed to verify the accuracy of the 
improved theory. However, these data are also needed to determine the accuracy 
with which correction factors derived from data at one Reynolds number can be 
used to predict pressures at another Reynolds number. This is particularly 
important for trailing edge control surfaces. 

A number of suggestions can be made for future wind tunnel tests in 
addition to those indicated above. Leading edge control surfaces should be 
tested; spoilers might also be considered. Very little data are available for 
these configurations. Models should be designed so that components can be 
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tested in their principal modes of motion. Complete models usually have 
moveable control surfaces but a moveable fin, horizontal tail, and engine 
pylon should also be considered to distinguish between component loads and 
interference loads. More oscillatory transonic data are needed. In two- 
dimensions, pitch data should be measured in addition to control surface data; 
in three-dimensions, data on both straight and swept wings are required. 
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CONCLUDING REMARKS 


CONCLUSIONS 

The basic conclusions arrived at as a result of the calculations and 
correlations presented are outlined as follows: 

( 1 ) One Set of Correction Factors Is Not Good Fo r All Modes 

Application of correction factors, determined from one mode, to other 
modes has not met with much success. Specifically, correction factors obtained 
using a pitch mode can not be applied to pressures due to control surface 
deflections. The converse is also true. In addition, application of correction 
factors, obtained using a pitch mode, to pressures due to a camber mode (and 
vice versa) have not proved to be very accurate either. Bergh and Zwaan, 
reference 6, on the other hand have concluded that correction factors calcu- 
lated using a pitch mode can be applied to a roll mode. These two modes, 
however, are similar. One of them has a constant angle-of-attack along the 
span while the other has a linearly varying angle-of-attack along the span. 

Further study is required to find out what types of correction factors 
are required for the various modes encountered in flutter and other dynamic 
aeroelastic analyses. The practical method of implementing such a correction 
procedure also requires further study. 

The fact that one set of correction factors can not be applied to all 
modes has certain implications for testing procedures. It may be that more 
rigid body types of modes will be required (e.g., pylon yaw, wing alone pitch, 
tail alone pitch, outer wing pitch, inner wing pitch, fuselage alone pitch, 
etc.) than are now considered. 

A second approach to the problem of obtaining one set of correction 
factors for both a pitch and a control surface mode was attempted using the 
"multiple mode" capability of the program. This capability allows the theory 
to be constrained to produce the correct lift and moment coefficient, etc., 
for each of several modes. The resulting span loading and/or pressures were 
not improved for either the pitch or control surface modes. 
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Even within a single mode problems can occur for different amplitudes. 

For instance high angle-of-attack flow fields can be basically different 
(separated) from those at low angles-of-attack. However correction factors 
obtained at low angles-of-attack can result in improved predictions for all 
angles-of-attack. The basic reason being that the viscous corrections, although 
smaller for the unseparated case, are still in the same direction as that for 
the separated case. 

(2) A Bending Moment Constraint Is Needed For Swept Mings 

Lift and pitching moment constraints are not enough for the swept wing 
case. A bending moment constraint is also required so that the loading is 
not shifted outboard to accommodate an aft shift in aerodynamic center of 
pressure. Without the constraint on the bending moment the correction factors 
will cause the wing loading to be moved toward the wing tip instead of moving 
the load aft along the chordline. 

(3) Correction Factors Can Be Extrapolated More Accurately To Other Mach 

Numbers Than To Other Frequencies 

When correction factors are determined at low Mach Numbers they are 
caused primarily by viscous effects. Since viscous effects exist at all Mach 
Numbers an increase in accuracy will result if the low Mach Number correction 
factors are applied to the high Mach Number cases. 

Extrapolation in frequency has not been as successful as extrapolation 
in Mach Number. For the two dimensional case studied it seems that extra- 
polation further than Ak r = 0.1 (based on the half chord and a Mach Number of 
0.85) will lead to a decrease in accuracy as k r is increased. It is believed 
that extrapolation in reduced frequency is more accurate at lower Mach Numbers. 

It appears that as the frequency is increased the viscous effects are 
reduced. This is an important fact and if steady wind tunnel results are to 
be used for correcting data then a good estimate of the reduction in the 
viscous effect must be known. One way to accomplish this without testing 
every configuration is to test a representative sample of configurations over 
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a range of frequencies and construct general trends to be used in conjunction 
with steady data to estimate the frequency effect on correction factors. 

(4) Qualitative Features Missing From The Theory Can Not Be Generated By 

Correction Factors 

Correction factors tend to produce quantitative changes to the theory and 
not qualitative ones. For instance, in transonic flow, the bulge in pressure 
at the shock location cannot be induced with correction factors if one did not 
exist in the basic theory. 

(5) Downwash Correction Factors Must Be Additive, Not Multiplicative 

It was found that postmulttplying correction factors did not work for 
control surface modes (they did work for pitch modes however). That is, scaling 
the downwash to reproduce the imposed experimental constraints (C^, C^, etc.) 
led to unusable results. Smoothing of the results was obtained by introducing 
correction factor modes; however, the levels of correction were still unreal- 
istic. An analysis was performed to see what downwash was required to produce 
the experimental pressures and it became evident that the downwash had to be 
corrected everywhere and not just on the control surface. This suggested an 
additive downwash correction. A new method was developed and executed 
successful ly. 

Downwash correction factors essentially reflect the physical fact that 
viscous effects tend to change the effective airfoil camber (and thus the 
downwash). The camber is changed over the entire airfoil. 

(6) Premulti plying And (New) Postmultiplying Correction Factors Are Equally 

Accurate 

In the cases studied the accuracy of the corrected theory is improved 
equally well (approximately) by either type of correction factor matrix. 

The downwash correcting factors (New Postmultipliers) are physically more mean- 
ingful if interpreted as viscous corrections while premultipliers are more 
meaningfully intepreted as compressbil ity corrections. 
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Since transonic experimental data reflect both viscous and compressibility 
effects a very accurate way to obtain correction factors, if data permit, is 
to combine pre- and postmultipliers together. First a postmultiplier (New) 
is developed at low Mach Number where viscous effects dominate. This correct- 
ion is then applied to the transonic case. The modified theory is then 
corrected further for transonic effects using a premultiplier. Correction 
factors produced in this way are more accurate than most. 

(7) New Transonic Method Useful But Requires Further Investigation 

Various methods of applying local Mach Number were tried. Simple proced- 
ures based on the substitution of the local steady Mach Number (or some average 
between surface and free stream) for the freestream value in the boundary 
conditions, kernel, and pressure equations have been tried. The results have 
only shown minor changes and have not even given qualitatively good results. 

A new method was developed at Douglas (under the McDonnell Douglas IRAD 
program) and is based on a transformation of the distance between sending and 
receiving points based on acoustic travel time between the two points. This 
method was implemented for the two-dimensional case and correlated in the 
present study. The results are encouraging since the predicted pressures are 
qualitatively similar to the experimental data. That is, the new method 
predicts a bump in the pressure which is centered at the shock wave location 
and predicts a lowering of the pressure forward of the shock wave. This bump, 
however, is forward of the experimentally observed bump, and is usually smaller 
in amplitude. 

The theoretically determined phase angles of the pressures are not in 
good agreement with the experimental data forward of the 40 percent point 
(for the case of a control surface rotation). One possible reason for this 
is the fact that the wave fronts emanating from points on the flap tend to 
move up and over the shock wave and arrive at the forward portions of the 
airfoil in nearly a horizontal configuration. On the other hand, in the theory, 
the paths of the wave fronts are assumed to be normal jto the free stream 

flow with the wave fronts vertical, and this difference causes the phase angles 
to be greater for the theory than for the data. Further investigation of this 
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discrepancy and its solution is required. It may be possible to use a more 
exact phase lag time computation in the theory, such as the one used by Tijdeman 
and Bergh in reference 11. 

When trying to decide which theory is best it is important to account for 
viscous effects. Without such a correction the new transonic method, designated 
as (Mj, seems best. However, when viscous effects are accounted for, the 
variation designated (M^) is best. 

With the current method of computation the new transonic method is 
probably not reliable past M^ = 0.90. This may not be too restrictive since 
most wings are swept which reduces the normal Mach Numbers to values lower 
than 0.90. 

(8) Overview of Conclusions 

The concept of correcting theoretical pressure or load distributions so 
that they reflect associated experimental data works well with the correction 
factor technique, especially if the proper experimental data are available 
(e.g. bending moments.) It was hoped however that a set of correction factors, 
once developed, would be applicable to a wide variety of other cases. The 
range of applicability however has not been as wide as hoped for. Success in 
extrapolating correction factors was obtained for Mach Number and to a limited 
extent for frequency. Attempts to apply correction factors to dissimilar mode 
shapes however has not met with much success. Therefore more than one set 
of correction factors is required. The use of several sets of correction factors 
to correct oscillatory aerodynamic generalized forces for use in dynamic 
aeroelastic analyses requires further investigation. 

Also concluded from the present analysis is that correction factors can 
not change the character of the load distribution. If a fundamental feature 
is missing from the theoretical loading then the correction factors will not 
make it appear. Thus theoretical methods must possess at least qualitative 
accuracy. 
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Recommendations for Further Studies 

Further studies may be profitably pursued in several areas. The successes 
and failures of the correction factor technique, presented here, furnish a 
guide to such studies. 

First, it seems advisable to exercise more of the various options in 
the present method and include more types of force data (integrated from pressure 
data) for some of the cases treated in this report. For instance such a case 
would be the Hertrich wing (refs. 14, 15). It would also be desirable to 
obtain new data similar to that obtained by Hertrich, in which both force and 
pressure data are available. It would be interesting to compare force data 
with integrated pressure data. 

Second, it is now clear that one set of correction factors is not suffi- 
cient for all deflection modes. Thus a method for including multiple sets of 
correction factors into the determination of generalized oscillatory aerodynamic 
forces for various modes is required. This method may require special testing 
procedures whereby each major component or subcomponent is systematically given 
a rigid body rotation. 

Third, the studies on viscous effects initiated in this report should 
be continued. Specifically the technique of determining theoretical camber 
lines that reproduce experimental pressure distributions seems valuable and 
could lead to a semi empirical method for viscous effects when combined with 
boundary layer theory. 

Fourth, the new transonic method illustrated in this report should be 
refined and extended. Initially the two-dimensional capability should be 
refined in the areas of; 1) pressure phase angle and, 2) unsteady shock wave 
motion. Subsequent to this a three-dimensional method should be developed. 

In addition, an investigation should be undertaken to explore the 
possibility of developing a semi empirical transonic method. The local steady 
Mach Number can be used as an adjustable parameter so as to produce the 
pressure distribution changes necessary to satisfy experimental constraints 

(lift moment etc.). 
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CORRECTION FACTOR COMPUTER PROGRAM 


Introduction 

As described in the Theoretical Development Section this method generates 
a set of correction factors that can be applied to a set of data (e.g., theo- 
retical pressure) such that the data satisfies certain imposed (e.g., experi- 
mental) constraints. 

For convenience this data will be referred to as pressure data since this 
is the most common application. However the correction factor procedure is 
not restricted to pressures and can be applied to other data sets (e.g., span 
loads, etc.). 

For this procedure it is assumed that one or more theoretical pressure 

distributions, AC p , (j = 1, number of pressure modes) are input. Associated 

with these pressures are: an area distribution, aA, a set of coordinates, 

(x, y, z), and a dihedral angle distribution y which are input via cards, 

tape or both. As an option the aerodynamic influence coefficient matrix, 

[A] = [D]" 1 , along with one or more normalwash distributions, w^, can be input 

in place of AC ; and as a matter of fact these must be input for postmulti- 
Pi 

plier correction matrices (i.e., correction factor matrices for the normal- 
wash). 

Constraint data (experimental data) are input as force or moment co- 
efficients. If a force coefficient, C g , is considered it is defined as 

C = - T aA aC t • f, (force coef.) (87) 

6 *\» u p a 


where c is a constant used to convert the dimensional sum into a coefficient 
form. For example if C 


C L then c is equal to the reference area 


The limits 


of the sum are also input to the program. The unit vector i g is in the 
direction of an input axis. A set of axes are input for use in the constraining 
and monitoring features of the program. Each axis can be input in one of two 
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ways; (1) a point and a direction or (2) by two points. The unit vector i- 
is calculated as follows: 

^ 

i. = i cosa + j cose + k cosy (88) 

a 

The unit vector "n is in the direction of the lifting pressure which is given 
in terms of the dihedral angle of the lifting surface. 

rf = ^(-siny) +"k(cosy) (89) 


where j and k are unit vectors in the y and z directions respectively and 
where a right handed system is employed where z is up, y is out the starboard 
wing and x is aft. 


If a moment coefficient, C g , is considered then it is defined as follows: 


C a = r I a A 
e 2f a 


("r x n) 


(moment coef.) 


(90) 


where 

r = (x - £^)i + (y - n^)j + (z - c^)k (91) 

and where are the coordinates of the first end point of the 

axis considered and where i is its direction. The constant c for the case 

d 

of C u has the dimensions of volume. 

The program has various other capabilities and one of these is its 
ability to monitor the corrected or uncorrected pressures. The integrations 
performed in equations (87) and (90) can be performed using data without reference 
to constraints. Thus if span loads are desired for data that has been corrected 
(or uncorrected) then the proper summations are activated in the program in 
a manner similar to that for constraining the data. 

The program also has the capability to use correction factor modes. That 
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is, the actual correction factors {e} are related to a set of modal coordinates 
{e} as follows: 


= [<f>] ( e> (92) 

The modal matrix, <f>, is either input directly by cards or certain built in 
modes can be activated. 

The program has the capability to limit the excursion of any or all 
correction factors. The upper and lower bounds are simply input for the corr- 
ection factors that are to be limited. If correction factor modes are used 
then the limits are placed on the modal coordinates, e, and not on the 
correction factors themselves. 

In addition to limits, a factor a is input for each constraint to 
indicate its "constraining power". The term a ranges from 0 to 1.0. If a 
is 1.0 the constraint, has full power and is 100 percent effective as a con- 
straint. If a is 0 then the constraining power is zero and the ^constraint" 
has no effect. For values of a anywhere in between the constraint is said 
to be an estimate. 

Finally, the program can be used to apply previously obtained correction 
factors to input pressure distributions. The program can also be used simply 
to monitor existing data without any constraints. 

One nomenclature problem which might cause confusion is the fact that 
the normalwash w is called W in the program, while the correction factors W 
are called CF. 
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Program Input 


The following table provides an overview of the card input data grouped 
according to their functions. in the program. The layout of the input sheets 
and a detailed description of each input item are also given following the 

table. 

Overview of Input Data 


ITEMS 

i 

CARD 

NO. 

WHEN 

NEEDED 

COMMENTS ' 

r 

CONTROL I 

DATA 

1,2, 
3, 4 

Always 

Header card, control dimensions 
and control flags. 

GEOMETRY 

AND 

PRESSURES 

5, 6 

If 

FLAGP=3 

Geometry data is input on one card 
per i, i=l, NP (NP=number of pressures). 
Pressures are input either 6 real 
numbers per card (when FLAGI=1), or 3 
complex numbers per card (when 
FLAGI=0) . Repeat pressure input per 
pressure mode, symmetric modes first, 
antisymmetric modes (if any) last. 

AXIS DATA 

7, 8 

Always 

Axis data, 2 cards per i, 

1=1 , NAXIS (NAXIS=number of input axes) 

CONSTRAINT 

DATA 

9, 10 
11, 12 

If 

NC t 0 

Constraint data is input in a minimum 
of four cards per i, i=l, NC 
(NC=number of constraints) 

MONITOR 

DATA 

13, 14 
15, 16 

If 

NM0N f 0 

Monitor data is input in a minimum 
of four cards per i , 

1-1, NM0N (NM0N= the number of monitored 
aerodynamic parameters) 

LIMITS 

DATA 

17 

If 

NELIMS^O 

Minimum and maximum limit values on e; 
repeat per i, 1=1, NELIMS 
(NELIMS=number of min. and max. limiting 
value pairs) 
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ITEMS 

CARD 

NO. 

WHEN 

NEEDED 

COMMENTS ‘ 

CORRECTION 

18, 19 

If 

Correction factor modes may be input 

FACTOR 

MODES 

20, 21 

NEM f 0 

according to two options depending on the 
flag TYPE (see detail description of data) 
either in cards 18, 19 and 21, or in cards 
18, 20 and 21. Repeat per i, i = 1, NEM 
(NEM=number of correction factor modes) 

DOWNWASH 

DATA 

22, 23 

If 

FLAGW=1 

Downwash data is input in a minimum of two 
cards per mode (see detail description of 
data). Input symmetric modes first, 
antisymmetric modes (if any) last. 


Computer program requires less than 200K OCTAL storage. 
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PUNCH IN ALL CAROS CASE PROGRAM NO. 
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PUNCH IN ALL CARDS CASE PROGRAM NO. 



DIRECTIONS FOR KEYPUNCH DO NOT PUNCH BLANK COLUMNS NO UNDERPUNCHES IN SIGN FIELDS. 




Description of Input Data 
Control Data 

These data items are required for all cases. They consist of a header, 
control numbers, flags, and tape (or scratch unit) numbers. 


CARD 

ITEM 

MNEMONIC 

~ DESCRIPTION 

1 

Header 

HEADER 

Alpha-numeric description of case in card 
columns 1 through 60 

2 

NP 

NP 

Number of AC p elements, where AC p may represent 
any type of quantity (NP £ 350) 

2 

NC 

NC 

Number of constraints to be applied to the 
AC p values (NC <_ 35) 

2 

NEM0DES 

NEM 

Number of correction factor modes if any (NEM <_ 100) 

2 

NELIMS 

NELIMS 

Number of input cards giving the minimum and 
maximum values of e (NELIMS <_ 100)- 

2 

NMONITOR 

NM0N 

Number of sets of monitoring data used to 
integrate AC p into aerodynamic parameters (NM0N _< 35] 

2 

NAXIS 

NAXIS 

Number of axes input for use in integrating the 

AC data into forces and moments for constraint 
P 

and monitoring purposes (NAXIS <_ 25) 

3 

FLAGB 

FLAGB, 

IFB 

FLAGB=0, correction matrix calculation 
FLAGB=1 , monitor data only 

FLAGB=2, apply input correction factor matrices 
to input pressure distribution 

3 

FLAGP 

1 

FLAGP, 

IFP 

FLAGP=0, geometry data and AC p are input from 
tapes; calculate premultiplying correction 

factors. (See Tape Description section for format) 

1 . 
FLAGP=1 , geometry data and D (inverse aero matrix, 

are input from tapes, W (normalwash) input either 

from tape or on cards (see FLAGW below), calculate 

post-multiplying correction factors 
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CARD ITEM 


3 FLAGT 


3 FLAGW 


3 FLAG I 


3 | IPRINT 


MNEMONIC DESCRIPTION 

FLAGP=2, input as for FLAGP=1 ; but calculate 
pre- multi plying correction factors 

FLAGP=3, geometry data and AC p input on cards; 
calculate pre- multi plying correction factors 

FLAGP=4, geometry data and D" 1 input from tapes, 
W input either from tape or on cards; calculate 
modified post-multiplying correction factors 

FLAGT FLAGT=0, weights for minimization process are 

absolute values of forces for unit deflections 

FLAGT=1 , weights are unity 

★ 

FLAGW FLAGW=0, normalwash matrix, W, is input from 

tape, if needed 

• • FLAGW=1 , normalwash matrix, W, is input on 

cards 

FLAGI FLAGI=0, AC p values are input as complex 

numbers (either from tape or on cards) 

FLAGI=1, AC p values are input as real numbers 
(i.e. not complex) 

IPRINT Detail print flag; 

IPRINT = 1, print rows of the 
SAI matrix, and rows of the 
SAN matrix (if any) 

IPRINT = 0, bypass printing 
of same 


* Cap. W is normalwash in the computer program where correction factors are 
called CF. 


73 



1 CARD 

ITEM 

MNEMONICS 

DESCRIPTION ^ 

am 

TS 

NMSYM 

Number of symmetric pressure modes 

m 

TA 

NMASYM 

(NMSYM, NMASYM < 10) 

Number of antisymmetric pressure modes 




Note that all data items in cards 2 through 

* 

4 are input as integers, right-justified 



i 

in their respective fields of ten card 
columns each (format 110) as shown on the 
input sheets. 


* Right justified means input ending in the last (or right-most) card 
column of the field. 
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Geometry and Pressure Data 


Data items defining the geometry of a case are usually available on 
tape; similarly, pressure data (if needed; see item FLAGP under Control 
Data) are usually input from tape. However, if this is not the case, these 
data items may be input from cards by specifying FLAGP=3, as shown below. 


CARD 

ITEM 

MNEMONICS 

DESCRIPTION 




The following two cards are input only if 
FLAGP=3. 

5 

x i 

X 


5 

y i 

Y 

x, y, z coordinates of pressure point i 

5 

z i 

Z 


5 

^i 

GMA 

Dihedral angle of pressure point i 

5 


DELA 

Area of box over which the pressure acts. 
Repeat card 5 for all points, i=l , NP 

6 

4C p 

DCP 

Array of the AC p values (lifting pressures) 
either 3 complex numbers per card (when 
FLAGI=0), or 6 real numbers per card (when 
FLAGI=1; see Control Data) 




The format used for cards 5 and 6 is 6F10.0. 
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Axis Data 


The following data items are required for all cases. These input data 
are used to describe an axis in space. Axes can be described by either two 
endpoints or by one endpoint and a set of direction cosines. These axes 
are used in the integration of the pressures into force or moment coefficients. 
Forces are resolved in the direction of the axis, while moments are taken 
about the axis. 


CARD 

ITEM 1 

MNEMONIC 

DESCRIPTION 

7 


IAX 

Axis number 


■ 



7 

Axis 

I FA 

IFA=0, axis endpoints are input 


type 


IFA=1, a point and direction cosines 
are input 

8 

Cl 

XII 



nl 

ETA1 

Axis endpoint coordinates 


Cl 

ZETA1 


] 

8 

£2/C0Sa 

XI2 

Second axis endpoint coordinates if IFA=0; 

I 

n2/cose 

ETA2 

direction cosines if IFA=1 


£ 2/cosy 

ZETA2 

1 



Card 7 format is 6110; 
card 8 format is 6F10.0 
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Constraint Data 


The correction factors modify the theoretical values of AC p by a 
minimum amount so that specified forces and moments are reproduced. For 
example, if the total lift is known experimentally, then several data items 
must be input specifying the actual value of the lift coefficient and 
describing the way the AC p values are to be integrated to obtain this co- 
efficient. The lift coefficient is then called a constraint on the 
theoretical data. 


CARD 

9 

ITEM 

Axis 

number 

MNEMONIC 

JAX 

— DESCRIPTION 

Number of the axis to be used for calculating 
the constraint force or moment 


F-M 

Flag 

IFF 

IFF=0, the constraint, C g , is a force 
in the direction of the axis; 

IFF=1, the constraint, C g , is a moment 
about the axis (right-hand-rule). 

Card format is 6110 

10 

6 

NDI 

6=1, symmetric pressure mode to be used; 
6=-l, antisynmetric pressure mode to be used 

10 

Press. 

mode 

MI 

Pressure mode number to be used with 
constraint C e 

10 

a 

AIT 

l 

Constraining effectiveness of C 0 ; 

0 < a < 1. 

If a = 1, C ■ is a constraint; 

Cl 

if a < 1, c e is only an estimate, and the 

resulting weighted (corrected) theory will 

only approximately reproduce C g • If a= 0» 

then C will not affect data, 
e 

10 

c 

CIT 

Constant used to nondimensional ize integrated 
data. 

If C fi is a force, c = Area; 

if C is a moment, c = Area x length, 

e 

J, — — — 
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CARD 

10 

ITEM 

C e 

MNEMONIC 

CIE 

DESCRIPTION" 7 

Experimental (or any other) constraint on 
the data. Card format is 2110, 4F10.0. 

11 

LIMI1 , 
LIMI2 

i 

LIMI(l) 
LIMI (2) 

Identification of a range of AC p values 
(or boxes) from LIMI 1 to L I MI 2 defining 
the limits of integration for the pressures. 
There may be as many sets of ranges input 
as needed. Card format is 6110. 

12 

. \ 

i 

-1 

-1 

The number -1 ; end indicator for the 
sets of data LIMI 1 , LIMI2 




Cards 9 through 12 are repeated for all 
constraints, i.e. NC times. 
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Monitor Data 


The following data are input only if the control data item NM0NIT0R 
has a value different from zero. These data are used for the integration 
of the AC values into some meaningful parameters as a check on the effect 
of the correction factors on theory, whenever FLAGB = 0. Since often it 
is desirable to monitor the unmodified data as well, the setting FLAGB - 1 
is designed to integrate the &Cp values into parameters without calculating 
the correction factors. Another monitoring option can be activated by the 
setting FLAGB = 2; in this case the correction factors are input from tape 
(FT16) saved in a previous run and the weighted ACp values are integrated 
into parameters as specified $>y the monitor data. 


CARD 1 

ITEM 1 

MNEMONIC 

________ 

13 

Axis 

number 

NAX 

Axis number used in the integration of the 

aC values into forces and moments 
P 

13 

F-M flag 

IFN * 

IFN=0, parameter to be determined is a 
force 

IFN=1 , parameter to be determined is a 
moment. 

Card format is 6110. 

14 

6 

NDN 

6=1, symmetric pressure mode used 
6=-l , antisymmetric pressure mode used 

14 

Press. 

mode 

MN 

Pressure mode number 

14 

'V 

a 

ANT 

Not used 

14 

% 

c 

CNT 

Constant used to nondimensional ize 
integrated data 

14 

LABEL 

LABEL 

Alphameric identifier of the integrated 
parameter (ten characters long) 

Card format is 2110, 2F10.0, 10A1 
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CARD 

ITEM 

MNEMONIC 

DESCRIPTION 

15 

LIMN1 , 
LIMN2 

LIMN(l) 
LIMN ( 2 ) 

Identification of a range of ACp values 
defining the limits of integration for 
the pressures. There may be as many sets 
of ranges input as needed. 

Card format is 6110. 

16 

-1 

-1 

The number -1; end indicator for the 
sets of data LIMN1 , LIMN2 



Cards 13 through 16 are repeated for 
all parameters, i.e., NM0NIT0R times 













Limits Data 


It is sometimes desirable to place a restriction on the range of values 
of e by specifying a minimum and a maximum bound on e. In this case the 
control data item NELIMS is input as the number of e limit pairs to be supplied, 
which are input as shown below. Note that this input (card 17) is omitted 
when NELIMS = 0. 


CARD 

ITEM 

MNEMONIC 

DESCRIPTION 

17 

LIMel 

LIMK(l) 


A range of boxes, or ACp elements, 


LIMe2 

LIMK(2) 


over which a limit is placed on t 

17 


EBMIN 

The minimum 



III 1 (I 



value allowed for e 


e 

EBMAX 

The maximum 



max 







Card format is 2110, 4F10.0 




Card 17 is repeated for all sets of 




ranges, i.e., 

- - _ 

NELIMS times. 
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Correction Factor Modes Data 

In many instances it is desirable to restrict the incremental correction 
factors {e} to a linear combination of a set of modes, { e } = [$] (e g > • The 
mode shapes [4>] can be input directly per box, and per mode, or the mode 
shapes may be selected from a set of functions. This set of data is input 
only if the control data item NEM is different from zero. 


W 


ITEM 


MNEMONIC 


DESCRIPTION 


18 

18 


18 

18 

18 


e MODE 
NUMBER 

TYPE 


n 

a 

b 


M0DEN0 

ITYPE 


NL 

AL 

BL 


Weight factor mode number 

TYPE=1 , use (x-a) n 

TYPE=2, use (y-a) n 

TYPE=3, use (z-a) n \ as 

TYPE=4, use exp[b(x-a) n ] ( mode 

TYPE=5, use exp [b(y-a) n ] 1 equation 

TYPE=6, use exp [b(z-a) n ] ) 

If TYPE = 0, the <|> values are input in 
card 19, and the following 3 items are 
not used 

Constants used in the mode equation 
Card format is 2110, 4F10.0 


19 


19 

19 

19 


*(J) 

J + 1 


PHI (J) 

J + 1 
PHI (J+l ) 


Card 19 is input only if TYPE = 0. 

Box (or element) number for which the 
<f>(J) applies 

The modal value of the J-th value of e 


) 


Another set of <|>-data 
Card format is 2(110, 2F10.0). 

Repeat as needed, 2 sets of data per 
card. Note that only the non-zero element] 
need be input. 







CARD 

ITEM 

MNEMONIC 

DESCRIPTION 




Omit card 20 if TYPE = 0. 

20 

IIML1 

LIML2 

LIML(l) 
L IML ( 2 ) 

Range of boxes or e's over which 
the current e-mode applies. There 
may be as many sets of ranges input 
as needed. 

Card format is 6110. 

21 

-1 

-1 

The number -1; end indicator of 
data set for the current e-mode 
(M0DEN0) 


Repeat cards 18 through 21 for all 
e modes, i.e., NEM times. 














Normalwash Data 


If the normalwash matrix [W] is needed (see control data FLAGP), 
and it is not available on tape, the control flag FLAGW must be input as 
1, and then the normalwash values are card input as shown below. 


CARD 

ITEM 

MNEMONIC 


DESCRIPTION 

22 

M0DE 

M0DE(J) 


Mode number for the current set of 
W values 

22 

6 

IDELW 


Symmetry flag to aid in identifying 
the mode; note that 6=1 type values 
are expected to precede the <5=-l type 
values 

22 

LIMW1 

LIMW2 

LIMW(l) 

LIMW(2) 


A range of boxes over which 
the W value applies 

22 

W 

WIN 


Normalwash, W, for the above range 
of boxes. 

Card format is 4110, 2F10.0. 

Repeat card 22 as needed. 

Note that only the non zero W values 
need be input. 

23 

-1 

-1 


The number -1; end indicator for 
the normalwash input data 
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Tape Description 

Program E1GC uses a minimum of four, and a maximum of twelve tapes 
and/or utility (scratch) units depending on the type of the case considered. 
In addition NPIT = 5 and NP0T = 6 are used throughout the program as the 
system input/output units respectively. These, as well as all tapes and 
utility units are defined in subroutine WEYT by means of a DATA statement 
specification under their respective names. The following table gives a 
summary of tape names and their use; the formats of those tapes that may 
be specified as input/output units are described in subsequent tables. 
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Summary of Tape Units 


NAME 

UNIT 

i WHEN 

NEEDED 

USER 

! SUBROUTINES 

DESCRIPTION OF CONTENTS 

NUTL1 

1 

Always 

WEYT, WSWA, 
SDBL, EPSJ 

Miscellaneous intermediate 

NUTL2 

■ 

Always 

WEYT, SDBL, 
DCPT, CEMN 

solutions 

NTSAIJ 

E 

If NC^O 

WEYT, SAIJ, 
DELC, SDBL 

SAI matrix rows 

NTSANJ 

fl 

If NM0N?*O 

WEYT, SAIJ, 
CEMN 

SAN matrix rows 

NTPHIJ 


If NEM^O 

WEYT, PHIJ , 
SDBL, EPSJ 

<j> matrix columns 

MASTSB 

B 

Always 

WEYT, SDBL, 
GINV 

S matrix rows 

NEWTSB 

10 

If NELIMS 
t 0 

WEYT, M0DF, 
GINV 

The modified S matrix rows 

NTGE0M 

11 

If FLAGP^3 

WEYT 

Geometry arrays; input tape 

NTDCP 

12 

Always 

WEYT, DCPB 

ACp matrix columns; 
either input tape or scratch 
unit depending on FLAGP 

NTAPW 

13 

| If FLAGP 

WEYT, WSWA 

W (normalwash) columns; 

input tape if NTAPW = 0, scratch 

unit otherwise 

NTAPDI 

14 

W 0, 3 

WEYT, DCPB, 
SDBL, DCPT 

Inverse downwash factor 
matrix [D]"^ 

NEWDCP 

15 

If FLAGI=1 
and FLAGP 

WEYT, DCPB 

Complex ACp columns, when 
ACp is input as a real matrix 

NTAPCF 

16 

- u 

If FLAGB^l 

WEYT, DCPT 

CF, the correction factor 
matrix; NTAPCF is output tape 
(or scratch unit) for FLAGB = 0 
cases; NTAPCF is an input tape 
for FLAGB = 2 cases 
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Input Tape NTGEgM 


RECORD 

WORD 

ITEM 

DESCRIPTION 

1 

1 

LENGTH 

Length of arrays in records 
2 through 6 (LENGTH = NP) 

2 

1 - NP 

X 

x-coordinate array 

3 

1 - NP 

Y 

y-coordinate array 

4 

1 - NP 

Z 

z-coordinate array 

5 

1 - NP 

GMA 

Dihedral angle (y) array 

6 

1 - NP 

DELA 

Array of box areas 


Input Tape (or Scratch Unit) NTDCP 


RECORD 

WORD 

ITEM 

DESCRIPTION 

1 

1 

• • NP 

Row dimension of the AC p matrix 
(column length) 


2 

NSYM 

Number of AC p columns for 
symmetric modes 


3 

NASYM 

Number of AC p columns for 
antisymmetric modes 

2 

1 - NP 

DCP 

aC column for first symmetric mode 
P 

• 

. 

• 

• 

• 

• 

• 

• 

• 

• 

• 

• 

1+NSYM 



AC p column for last antisymmetric 
mode* 

+NASYM 

. _ 1 

1 - NP 

i 

DCP 


* Note that if NASYM = 0, the last &C p column refers to last symmetric mode 
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Input Tape (or Scratch Unit) NTAPW 


RECORD 

WORD 

ITEM 

DESCRIPTION 

1 



1 

NP 

Row dimension of the W 
(normalwash) matrix 


2 

NSYM 

Numbers of W columns for 
symmetric modes 


3 

NASYM 

Number of W columns for 
antisymmetric modes 

2 

■a 

W 

W column for first symmetric 
mode 

• 

• 

• 

• 

• 

• 

• 

* 

• 

• 

• 

• 

1 + NSYM 
+NASYM 

1 - NP 

W 

W column for last anti- 
symmetric mode (if any) 
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[ WORD ITEM 

DESCRIPTION 

1 

NP 

Row dimension of matrix DI 

2 

NP 

Column dimension of matrix DI 

3 

NC0L2 

NC0L2 = NP if both symmetric 
and antisymmetric DI matrices are 
on tape; NC0L2 = 0 otherwise 

1 - NP 

DI 

First row of DI, the inverse 
downwash factor matrix, [D 
for synmetry 



The following records may be 
omitted when antisymmetric modes 
are not desired. k 


1 _ np DI First Di-row for antisymmetry 


1 - NP 


DI 


Last Di-row for antisymmetry 














Input/Output Tape NTAPCF 


RECORD 

WORD 

ITEM ] 

DESCRIPTION | 

1 

1 

C0DE 

Alphameric identifier of tape, 

4 characters in length, left justified; 
C0DE = PRE, for pre-multiplier 
cases, C0DE = P0ST for post- 
multiplier cases 

2 

1 

LENGTH 

Length of array CF. LENGTH 
should be equal to NP. 


2 

NMSYM 

Number of symmetric modes for case 

1 

3 

NMASYM 

Number of antisymmetric modes for case 

Note that the last two items are not 
used when tape NTAPCF is an input tape 

3 

HD 

CF 

Array of the complex correction factors 
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Test Cases . - The use of the program will be illustrated by two test 
cases. The first will be a premultiplier and will exercise most features 
of the program so that their use can be illustrated. The second test case 
will illustrate the use of the rvew postmultiplier, tape input, and downwash 
input on cards. 

The theoretical pressures are taken from a two-dimensional analysis of 
an airfoil with a 25 % chord flap. The new transonic procedure discussed 
previously will be used for the airfoil operating at a Mach Number of 0.875 
and a reduced frequency of 0.0. Figure 47 illustrates the geometry, pressures 
and axes data for the airfoil for control surface rotation (Mode 1) and 
pitch (Mode 2). Also shown on the figure are the theoretical and experimental 
values of c £ , c ml/4 and c h]/4 for mode 1 and c £ for mode 2. The experimental 
values are used as the constraints. An experimental value for c £ for mode 2 
is not available thus an estimate is given in its place in the figure. 

Test Case 1 . - Table III presents the input cards for the first test case. 

The number of pressures, NP, is 19; the number of constraints, NC, is 4. For 
this case 19 correction factor modes, <j>, will be used, thus NEM = 19. In 
addition limits will be placed on the values of e. These limits will be 
described by one card thus NELIMS =1. The number of axes, NAXIS, is 3. 

The program is able to monitor the corrected data, and in this test case the 
number of coefficients to be monitored, NMON, is 4 and they are c^, c ml ^ 4 and 

c h3/4 f0r m0de 1 and C * f ° r m0de 2 ‘ ThUS thC monitored coefficients should 
reproduce the input constraints. This, in fact, is the case as the output 

shows in Table IV. 

Since correction factors are to be calculated rather than data monitored 
only, FLA6B = 0. Also since the geometrical data and pressure data are to be 
card input and a premultiplier is to be calculated FLAGP = 3. The usual 
weight factor T, (the absolute value of the force on an element) is not used, thus 
FLAGT = 1. Normalwash values are not input thus FLAGW = 0. Only real values 
of pressure are used thus FLAGI = 1 ; the detail print flag is input as 
IPRINT =1. In this example there are two modes (call them symmetric) thus 
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NMSYM = 2 and NMASYM = 0. This marks the end of the control data. 

The geometry data are taken from figure 47 and are given on cards 
designated as type 5. The 1/4-chord point of each box is input along with 
its area, aA = aX. The pressures at each 1/4-chord point of each box and for 
each mode are taken from figure 47 and are given on cards designated as type 6. 

The axis data are encountered next. IAX identifies the axis number and 

1 FA identifies how it is input (whether by two points or a point and a dir- 
ection). In this case a point n^, and direction (cosa, cose, 

cosy) are input thus IFA = 1. These points and directions are taken from 
figure 47 and are input on card designated as type 8. 

The constraint data is next. Input are four constraints c £ , and 

c L , (J for mode 1 and c 0 for mode 2 taken from the experimental values of these 
parameters given on figure 47. Each constraint has a 9 and 10 type card. 

JAX identifies the axis to be used with the constraint (axis 1 for c £ , axis 

2 for c ml/4 and axis 3 for c h3/4 ). The flag IFA identifies the coefficient 
type tote calculated whether the force type (IFA = 0) or moment type (IFA = 1). 
The terms MI and NDI denote the mode to be used. In this case modes 1 and 2 
are symmetric. The constraining power AIT is taken as 1.0 for the constraints 
of mode 1 to ensure a full constraint. However AIT for c £ of mode 2 is taken 
as .95 since this is an estimate. The nondimensionalizing constant CIT is 

the chord for the c £ constraint and the chord squared for c ml/ , 4 and c h3 / 4 . 

The limits of integration LIMIT . LIMI2 span the entire surface for c £ and 

c , ,, (from box 1 to box 19) but only range over the control surface (box 13 
ml/4 v 

to box 19) for c h3 ^ 4 - 

The monitor data found on card types 13, 14, 15, 16 are almost identical 
to that of the constraint data because in this case c £ , c m ^ 4 and c h3 ^ 4 are 
the parameters to be monitored. Of course they could be any quantity or 
for that matter no quantities if monitoring is not desired. The only real 
difference between monitor data and constraint data is that an alpha-numeric 
identifier is input in place of the constraints for the monitor data. 
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As an example of the use of limiting values on e, card type 17 is input 
for this test case. Specifically it is required that 


-0.7 <_ e 1.5 

hold for all values of e, 1 through 19 {LIMK1 = 1, LIMK2 = 19). 

As a simple example of the use of correction factor mode shapes, <j>, an 
identity matrix will be used; 

[♦] = ru (itype = 0) 

Card types 19 and 21 are used to input these modes. 

The program output for this case is given in Table IV. The printed out 
put, which fits on 8 1/2 x 11 sheets, contains most of the input. Integration 
matrices are then printed along with other intermediate steps in ithe process 
of solution. At the end of the printout a summary of the geometry data, 
incremental correction factors, e, and modified pressures are printed. Next 
are the correction factors e + 1 and finally the aerodynamic parameters, 
calculated using the modified pressures, that have been monitored by the 
program. 

Test Case 2 . - Table V presents the input sheets for the second test case. 

This test case is the same as the test case 1 with the following exceptions: 
(1) the geometry and [D]" 1 are input from tape; (2) a new postmultiplier is 
developed; (3) one mode is used with three constraints; (4) correction factor 
modes are not used and (5) limits on incremental correction factors, e, are 
not imposed. 

For this case changes from test case 1 occur in the control data (cards 1 
through 4). First, no correction factor modes (NEM = 0) are to be used. 
Second, the card giving limits on e is omitted, thus NELIMS = 0. Third, [D] 
and the geometry are to be read in on tapes and the new postmultiplying 


93 


correction factor is desired, thus FLAGP = 4. In this case normalwash values 
are to be card read and so FLAGW =1. Also only one mode is to be used 
(control surface rotation), thus NMSYM =1. The geometry data remains the 
same as in test case 1. 

Finally the normalwash is input on card type 22 . The mode is a control 
surface rotation, thus W =1.0 over boxes 13 through 19, i.e., LIMW1 - 13, 
LIMW2 = 19. The program output is given in Table VI. 


* Remember W is downwash in this computer program 
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TABLE III 
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CONTROL FLAGS 


TABLE IV 


0UTPUT LISTING - TEST CASE 1 



i 

i 




GO 




a 




CL 


1 


< 


1 


K 




►— 





^ sO 

GO 


CL 


Z 


h- 


o 


3 

II It II II M 



C 


GO 



UJLULLOi 

z 


K 

a cl cl a i/> 

LU 


D 

<1 <r ^ 

X 

QN ^ 0- p- 4 m 

CL 



rH 

Z 

K 

C 


►—1 

LUO 




>a_ to< 

_J 

II M M It II II 

u. 

oco oau 

c 

GO 

o 

K 1 ll 

OC 

X oo 


uj< > * 



K 

XI- ZCX 

Z 

uox 

GO 

O-j •—a' 

o 

ClOLL'LuX< 

»— t 

LULU 1 O 

o 

ZZZZZ2 

-J 

Oqjoo 


k 


98 



CARD-PEAD OPTION IS SPECIFIED 
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SYMMFTRIC DELTA-CP MATRIX 
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THE 3 SETS OF INPUT DATA FUR ALL AXES 
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THE 4 SETS OF INPUT DATA FOR Mf N I TOR I NG 

AXIS NO ( N ) FLAGF(N) DELTA M L I Ml LIM 2 A-WIG(NI C-WIG(N) LABEL 
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S A I MATRIX ROWS 
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CORRECTION FACTORS ** 
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SYMMETRIC CF-N MODE 
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TABLE V 

INPUT DATA - TEST CASE 2 
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19-1*0 NORMALWASH DATA 


TEST CASE NH. 2 NEW POST-MULT., T APE INPUT, M-.875 


TABLE VI 


0UTPUT LISTING - TEST CASE 2 
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THE 3 SETS OF INPUT DATA C DR ALL CONSTRAINTS 

AXISNO(I) FLAGF(I) DELTA M L I M 1 LTM 2 A— W I G C I ) t-WlGU) C-FCI) 
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AI MATRIX POWS 
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IS CARn INPUT 
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THF l COLUMNS n«= TMP OF! TA-OP-RAP MATRIX 
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Subroutine Description 


The computer program for generating correction factors (E1GC) consists 
of twenty subroutines. The MAIN of this program reads and writes the header 
card and reads the control dimensions for a case; the latter are used for 
dimensioning most of the complex arrays that are passed into Subroutine 
WEYT via the argument list. Subroutine WEYT is the actual working main 
of the program, which calls all the major subroutines, supplying these with 
the necessary information via their argument lists. The following is a 
detailed description of all subroutines of program E1GC including their 
flow charts, where applicable, given In alphabetical order. The computer 
program is written in the FORTRAN IV programming language. 


142 



SUBR0UTINE CEMN(NP0T, IGQ, MQDE. NTAPSA, NP, NMON, LABEL, NUTL, 
SAI , DCPTIL, CE) 


Functional Description 

This subroutine integrates the corrected pressures, AC p , into 

coefficients, C g , which are used to monitor the results (Se§ Eq. (3)). 

The integration procedure is identical to that required for obtaining 

the imposed constraints. 

(Cl = [S] (AC } 

e K e 

The coefficients C g are part of the printed output. 



Description of Argument List 


NP0T 

IG0 

M0DE 

NTAFSA 

NP 

MM0N 

LABEL 

NUTL 

SAI 

DCPTIL 

CE 

Calling 


Data set number of the system output data set 

1 for symmetric modes 

2 for antisymmetric modes 
Mode number 

Data set (tape) number of tape containing the 
integration matrix [S.] 

Number of rows in the AC p matrix 

Number of integration rows used for monitoring 

Alphanumeric label describing the aerodynamic 

parameters 

Data set (tape) number of tape containing columns 

of the weighted pressures, AC, p 

A row of the integration matriS [S] 

A column of the weighted pressures {AC } 

^e 

A column of the aerodynamic parameters (C } 


Subroutine WEYT 
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SUBR0UTINE DCPB ( NTDCP , NTAPW, NTAPDI, 160, IFP, IFW, NR0W, NC0L, 
NMAX. DCP, C0L, W0RK) 


Functional Description 

This subroutine computes the theoretical pressure distribution if 
it is not input. Specifically 


{AC } = [D]" 1 (W) 


,-l 


where W is the normal wash and [D]~ 1 is the inverse of the aerodynamic 
influence coefficient matrix. This corresponds to Equation (1) where [D]' 
= [A]. AC is called DELCPB in this subroutine. 
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Description of Argument List 


NTDCP 

Tape number containing the matrix of pressure 
coefficients, [aC ], in column order 

NTAPW 

t 

Tape number containing the normalwash matrix, 
[W], in column order 

NTAPDI 

Tape number containing the inverse-D matrix, ([A]) 
[DI], in row order 

IG0 

1 for symmetric modes, 

2 for antisymmetric modes 

IFP 

Control flag (see input flag FLAGP). 

IFP = 0, 2, 3 means premultiplying correction 
factors, IFP = 1, 4 means post-multiplying 
correction factors 

IFW 

Normalwash flag. IFW = 0 means normalwash 
is tape input (if any), 

IFW = 1 means normalwash is card input 

NR0W 

Number of rows in the AC matrix 

P+ 

NC0L 

V 

Number of columns in the ACp matrix 

NMAX 

t 

Maximum number of columns in the aC matrix 

DCP 

Z 

One column of the AC matrix (complex) 

rf 

C0L 

w 

Temporary work array (complex) 

W0RK 

The NROWxNMAX complex array containing the aC matri 

v t 


Calling Subroutine 


WEYT 






SUBROUTINE DCPT(NP0T, FLAGB , IG0, M0DE, NP, NSCRCH, NUTL, NTAPDI ^ 
NTAPW, NTAPCF. X. Y. Z, GMA, D ELA, NMAX, NEM, W t . 
PI T EPS. DCPBAR . DCPTIL, W0RK, EB) 


Functional Description 

This subroutine modifies the theory with the calculated correction 

factors. If a premultiplier is used the theoretical pressure, AC p ^, is 

modified to produce the modified pressures AC (see eqs.(2) and (5)). 

r e 

{ 4c Pe } = 

If a postmultiplier is used then the downwash, W, is modified to produce 
the corrected pressures AC p (see eqs. (5) and (28)). 

|AC p | = [D ]’ 1 H+eJ {W} 

If the new postmultiplier is used then 

| AC p | = [D ]" 1 (W + [♦] (el) 

[♦] = [♦] r*j 

U> = [<t>] T ( aC } 

p t 

(See eq. 65) 

Also the correction factors, CF, are written on tape where 
CF = 1 + e 
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Description of Arguments 


NP0T 

FLAGB 

IG0 

M0DE 

NP 

NSCRCH 

NUTL 

NTAPDI 

NTAPW 

NTAPCF 

X 

Y 

Z 

GMA 

DELA 

NMAX 

NEM 


Data set number of the system output data set 
Option flag for correction matrix calculation and/or 
monitoring of data 

1 for symmetric modes 

2 for antisymmetric modes 
Mode number 

Number of row elements in the ACp matrix 

Data set (tape) number containing the <f> matrix in 

row order (for FLAGP=4 cases only) 

Data set (tape) number on which the AC columns 

p e 

are saved 

Data set (tape) number containing the D"^ matrix rows 
(if needed) (D ^ = A) 

Data set (tape) number containing the W matrix columns 
(if needed) 

Data set (tape) number on which the matrix of correction 
factors, CF, is saved in column order 
x coordi antes \ 

y coordinates ( of the pressure points or the AC 

( h 

z coordinates 1 

Dihedral angle array of the boxes over which the 
pressures act 
Array of box areas 

Column dimension of the two-dimensional complex 
array W0RK 

Number of correction factor modes 
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w 

A column of the W matrix (complex) 

DI 

A row of the D' 1 matrix (complex), (£A] matrix) 

EPS 

e array 

DCPBAR 

A column of the aC matrix (complex) 

Pf 

DCPTIL 

A column of the ACp matrix (complex) 

W0RK 

G 

Two dimensional complex array containing the aC 

EB 

e array (e = <pe) 


matri 

t 


Calling Subroutine 


WEYT 
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ow Chart 
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SUBROUTINE DELC(NTAPE, NP0T, NC, NP. NM0DE , NMAX, CIE, PCI, SAI, W0RK) 


Functional Description 

This subroutine forms the difference between the theoretical, C^., and 
the experimental (constrained) C g coefficients (see Equations (9) and (10)). 

{aC 0 } = {C e } - [S] UC p } 

t 

It also prints out C t ( = [S] (AC } ) and AC. 

t p t 
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Description of Argument List 


NTAPE 

NP0T 

NC 

NP 

NM0DE 

NMAX 

CIE 

DC I 
SAI 
WORK 


Tape containing rows of the integration matrix, [S] 
Data set number of the system output data set 
Number of constraints applied to AC p values 
Number* of row elements in the AC p matrix 
Number of modes 

Maximum number of columns in the two-dimensional 
W0RK array 

Array containing the input values 

C £ (experimental constraints) 

The [AC e ] matrix 

A row of the integration matrix [S.J 

The NP by NMAX complex array containing the AC matrix 

H t 


Calling Subroutine 


WEYT 







Functional Description 


This subroutine compares the correction factors, e, with the input 

limits e e . If any e falls outside of the limits it replaces e 
min’ max J 

with the closest limit. (The values of e are correction factors if 
correction factor modes do not exist). This subroutine forms the final 
correction factor array e (see paragraph below Eq. (60))* 
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Description of Argument List 


NP0T 

NELIMS 

NP 

NS 

LIMK 


JARR 

NSM0D 

EBMIN 

EBMAX 

EB 

ELIM 


Cal 1 i 


Data set number of the system output data set 
Number of input cards for EBMIN and EBMAX - see below 
Number of row elements in the ACp matrix 
NS = NP + NC when NEM < NP 
NS = NEM + NC when NEM > NP 

A two-dimensional array containing the first- and last box 
numbers that define a range of boxes (or AC p ) over which 
a limit is to be placed on e 

Array of the box numbers for which the e values are 
modified 

The number of e values which are modified due to the 
limits placed on these 

The minimum- and maximum value allowed for the values 
of e for boxes (or AC p ) in the range defined by LIMK 
Array of the calculated e values 
Array of the e values that were modified due to the 

Sin- Sax restrictions 
Subroutine WEYT 
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Flow Chart 
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SUBR0UTINE EPSJ ( NTPH I J . NP, NEM. NS, EB, EPS, PHI) 


Functional Description 

This subroutine relates e to e, as in Equation (53). 


(el = [<i>] (e) 

where [4>] are correction factor modes and where 

e = t £p for premultiplying correction 

( e for postmul iplying correction 
w 


factors 

factors 
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Description of Argument List 


NTPHIJ 

Tape number containing the ♦ matrix 

NP 

Number of row elements in the 4» matrix 

NEM 

Number of correction factor modes 

NS 

NS = the greater of (NP+NC) and (NEM+NC) 

EB 

Array of the e values 

EPS 

The final e array 

PHI 

A column of the matrix of weight factor 
shapes, <|> 


Calling Subroutine 


WEYT 










SUBR0UTINE GINV(NP0T, NTAPSB, NC, NS, NX, DC, EB, B, S, SBB) 


Functional Description 

This subroutine provides a general inverse of the following set of 
equations: 

i - 'V A 

NC J {AC e l = [S]{e)jNS 
When NC = NS (Direct Solution) 

* * _i 

{ e } = [S] UC e ) 

When NC < NS (Minimization Solution, £e = min (see Eq. (20)) 

{£> = [S] H {X} 

{X} = CCS] [ s ] H ] _1 {AC e > 

2 

When NC > NS (Least Squares Solution £ AC fi = min.) 

(e) + [[S] H [S]] 1 (AC e ) 

{X} = [S] H {AC e > 

In the above; 

cs] - [!] isr^ ]"’ 

where [S] given in Eq. (51). The term e is given in Eq. (22). The super- 
script H indicates the Hermitian transpose. 
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Description of Argument List 


NP0T 

NTAPSB 

NC 

NS 

NX 

DC 

EB 

B 

S 

SBB 


Data set number of the system output data set 
Tape number containing the S matrix 

* 

Number of constraints; number of rows in the S 
matrix (SBB) 

NS = the, greater of (NP+NC) and (NEM+NC); 

fV 

number of columns in the S matrix 
NX=NS if NEM=0, NX=NEM + NC otherwise 

The complex AC g array 

The complex array {e}» output of subroutine GINV 
An array of intermediate solutions (complex) 

A complex two-dimensional work array of dimension 
NC by NC 

% 

The complex NC by NS matrix, S 


Calling Subroutine WEYT 

Called Subroutine MIS2, the standard IBM system subroutine 

for solving complex matrix equations 
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SUBR0UTINE MATM(NT, 160, NR, NC. NMAX, A, C, B 


Functional Description 

Subroutine MATM is essentially a matrix multiplication routine. 

It obtains the DI matrix ([A]) rows from tape NT and the W matrix from 
the two-dimensional array B. The results of the matrix multiplication, 
AC p , are saved in array B which is returned to the calling routine, 
dcpS, via the argument list. 

Description of Argument List 

NT Tape number containing the inverse-D matrix 

IG0 1 for symmetric modes 

2 for antisymmetric modes 

NR Number of rows in the AC p matrix 

NC Number of columns in the AC p matrix 

NMAX Maximum number of columns In* the AC matrix 

*t 

A A row of the inverse-D matrix, [a] 

C Complex work arry 

B Two-dimensional complex array in which the AC p 

matrix is stored t 


Calling Subroutine 


168 



SUBROUTINE M0DF(NC, NS, MASTSB, NEWTSB, JARR, SQRTT , ELIM, 

SBB, PCI, DCM0D) 

Functional Description 

When some of the values of l have exceeded their limits and have 
been replaced by the limit values, these new values of e (called e d in 
Equation (56) are then considered fixed and known. However the constraints 
are now not satisfied and a change in the constraint AC e , i.e., AC mod , 
is calculated (see Equation (59)). 

AC mod = AC e " {£ d } 

Since the new values of e, i.e., e d , can not influence solution further 
the S matrix must be changed to delete the influence of e d . Thus the 
elements of S that give the influence of e d , i.e., S d , must be eliminated 
resulting in § u . This subroutine forms § u , or in the notation of the computer 
program [S^]. 
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Description of Argument List 


NC 

NS 

MASTS B 

NEWTSB 

JARR 

SQRTT 
ELIM 
SBB 
DC I 
DCM0D 


Calling 


Number of constraints - dimension of the complex 

arrays DC I and DCM0D 

Dimension of the complex array SBB 

Tape number containing SBB arrays (i.e., rows of the 

S matrix) 

Tape number containing the modified SBB arrays (i.e. 

rows of the matrix) S TOd = S u 

Array of the element numbers for which the S-values 

are replaced by zeroes 

/f7 - see Equations (23) and (34) 

J 

e, . array of the modified e values 
limj 

Complex array containing rows of the S matrix 
Complex array containing AC g 
Complex array containing aC^ 


Subroutine WEYT 
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SUBR0UTINE PHIJ(NPIT, NP0T , NTPHIJ, 



Functional Description 


This subroutine forms the correction factor modes. If <|i is input 
element by element (TYPE = 0) then this subroutine simply arranges the 
data into arrays. If TYPE = 0 then modes are calculated as follows: 


< x j " a /* TYPE = 1 

n. 

(y d - V 4 2 

n. 

(2 j - a *> 3 

exp[b^(xj - a^) 4 ] 4 

exp[b^(yj - a £ ) £ ] 5 

exp[b £ (Zj - a & ) A ] 6 


where a b n are input per mode and where <{»... = 0 over boxes are not considered. 


172 


Description of Argument List 


NPIT 

NP0T 

NTPHIJ 

NEM 

NP 

K0DE 

M0DES 

X 

Y 

Z 

PHI 


Data set number of the system input data set 
Data set number of the system output data set 
Tape number containing columns of the NP by NEM 
<j> matrix 

Row dimension of the <j» matrix 
Column dimension of the <f> matrix 
-1 (end indicator of card input sets) 
not used 

i 

y ' coordinates of the pressure points or 

z 1 of the ACp's 

Complex array containing one column of the <t> matri 


Calling Subroutine 


WEYT 













SUBR0UTINE P0SN(NT, IG0) 


Functional Description 

This subroutine positions tapes of a certain uniform format for 
reading; see Tape Description for NTDCP, NTAPW and NTAPDI. 

Description of Argument List 

NT Tape number to be positioned for reading 

IG0 1 for symmetric modes, 2 for antisymmetric modes 


Calling Subroutines 


DCPB, DCPT, MATM, SBAR 


SUBR0UTINE RECD(NTAPE, A, N) 


Functional Description 

This subroutine reads arrays pf real numbers A of length N from tape 
NTAPE one record at a time. It is used for the reading of the geometry 
arrays when these are input from tape NTGE0M. 

Description of Argument List 

NTAPE Tape number 

A Array to be read from tape 

N Length of array A 


Calling Subroutine 


WEYT 


SUBR0UTINE SAIO ( NPIT ■» NP0T, NTSAIJ, NTSANJ, NC, NP, NM0N, NAX IS, 
AIT t CIE. X. Y, 1 , CG, SG, DELA, FLAGA , FLAGF, 


K0DE, IPRINT, LABEL, SAI) 

Functional Description 

This subroutine sets up proper argument lists for SROW so that 
integration matrices> [S], can be calculated for both constraining 
and monitoring purposes. 



Description of Argument List 


NPIT 

NP0T 

NTSAIJ 

NTSANO 

NC 

NP 

NM0N 

NAXIS 

AIT 

CIE 

X, Y, Z 
CG, SG 
DELA 
FLAGA 


FLAGF 


Data set number of the system input data set 

Data set number of the system output data set 

Tape number containing the integration matrix rows, 

SA. for constraints 
1 J 

Tape number containing the inegration matrix rows, 

SAnj, for monitoring 

Number of constraints 

Number of AC p elements 

Number of sets of the monitoring data 

Number of axes 

Constraining effectiveness of the experimental data. 
Experimental (or any other) constraint on the data, C 
Coordinates of the pressure points (boxes) 

Cosine-, sine of box dihedral angles 
Box areas 
Axis flag 

= 0, axis endpoints are input 
= 1, direction cosines are input 
Force/moment flag 

. o, c Is a force in the direction of axis 

_ i , c is a moment about axis 
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K0DE = -1 

IPRINT Detail print flag 

I PRINT = 1, print SA . • 

* J 

and SA . rows 
nj 

IPRINT = 0, bypass print 

LABEL Alphameric identifier of the integrated parameters 

SAI Complex array containing a row of either one of the 

integration matrices SA^. or SA n j 


Calling Subroutine WEYT 

Called Subroutine SR0W 


179 



Flow Chart 
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SUBR0UTINE SBAR(NTSBIJ, NTAPDI , NC, NP, NS, FLAGP, FLAGT , FLAGW, 
I. I G0 , SQRTT, AIW, SAI, PI, W, DELCPB, SBI ) 


Functional Description 


This subroutine solves for the matrix SBI where 

-1 

sbi = [s] r Jf* j 

The matrix [S] contains all the capability of the program except modes 
and limits. This capability is outlined in Eqs. (45), (26) and (9) for pre- 
multipliers and (45), (39) and (30) for postmultipliers. The weights T are 
defined below Equation (47). 


SBI 


I 

I 


SA.. when j = 1, 2***NP, i = 1, 2***NC 

* J 


WT^ when j = NP + i 
0 otherwise 

where NP = number of pressure values and NC = number of constraints. 


SA ij 


WT i 


t j 


when ACp available 


/ s ftjj iCp .y /tt 

! y SA* DI 0 W. / /T- when AC not available 

/ *» ^ •> / 3 pj 

( 


% % f\j 

(1 - a-)/a. a. > 0.0001 


i " i 


m 

V 


10" 




a. <_ .0001 


( 

l 


| AC | AA FLAGP \ 1 


or 


1 , 0 


FLAGT = 0 


FLAGP = 1 
FLAGT = 1 
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Description of Argument List 


NTSBIJ 

NTAPDI 

NC 

NP 

NS 

FLAGP 

FLAGT 

FLAGW 

I 

IG0 

SQRTT 

AIW 

SAI 

DI 

W 

DELCPB 

SBI 


Tape number containing rows of the integration matrix 
without weight factor modes, S 

Tape number containing rows of the inverse-D matrix 
Number of constraints 
Number of AC p elements 
Length of the SBI matrix rows 

Option flag for &C p and/or pre- or post-multiplying 

correction factors 

Option flag for weights 

Option flag for normalwash input 

An intermediate index 

1 for symmetric modes 

2 for antisymmetric modes 
= / T ; see equations 
Constraining effectiveness a. 

A row of the integration matrix S'. 

• J 

A row of the inverse-D matrix 
A column of the normalwash matrix 

A column of the aC matrix (lifting pressure coefficients, 
p e _1 

either input, AC p , or computed as [D] (W1 ) 

A row of the integration matrix without weight factor 

, r?n r /t* i“1 t+ fT for constraints 
modes, [S] l/T* J , T* = I ^ 

\ for estimates 

/ a 


Calling Subroutine SDBL 
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Flow Chart 








SUBROUTINE SDBL(NSCRCH, NUTL, MASTSB, NTPHIJ, NTAPW, NTSAIJ, 

NTAPDI . IG0, FLAGW, FLAGP, FLAGT , NC, NP, NS, 

NEM. SQRTT , AIT, DELA. SBB, SBI, SAI, PI , W, PHI, DELCPl 

Functional Description 

B 

This subroutine calculates [S] described in Equation (54). The quantity 
calculated in this routine includes estimates and thus 


[S] = [SBB] 


[S] 


<t> 1 0 

- -I - 


0 »I 
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Description of Argument List 


NSCRCH 

Tape number containing columns of the <p matrix 


(if any) 

NUTL 

Utility (scratch) tape number 

MASTSB 

Tape number containing the S matrix rows 

NTPHIJ 

Tape number containing the $ matrix columns 

NTAPW 

Tape number containing columns of the normalwash matrix 

NTSAIJ 

Tape number containing rows of the integration matrix, SA 

NTAPDI 

Tape number containing rows of the inverse-D matrix 

IG0 

1 for symmetric modes 


2 for antisymmetric modes 

FLAGW 

Option flag for normalwash input 

FLAGP 

Option flag for AC p input and/or pre- or post-multiplying 


corrections 

FLAGT 

Option flag for weights 

NC 

Number of constraints 

NP 

Number of AC p elements 

NS 

= max (NP+NC, NEM+NC) 

NEM 

Number of correction factor modes 

SQRTT 

/ T. , see equations 

J - 

AIT 

Constraining effectiveness a. 

DELA 

Box areas 

SBB 

A row of the S matrix (integration matrix with weight 


factor modes) 

SBI 

A row of the [S] T v/t*J~^ matrix 

SAI 

A row of the S... elements 
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DI 


A row of the inverse-D matrix 

A column of the normalwash matrix 

A column of the <|> matrix (weight factor mode shapes) 


W 

PHI 

DELCPB A column of the AC p 


Calling Subroutine WEYT 

Called Subroutine SBAR 


matrix 

t 
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Flow Chart 
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SUBR0UTINE SR0W(FLAGA. FLAGF, XII, ETA1 , ZETA1 , CG , SG, CTIL, 
X. Y. Z, DELA, LIMI, UMAX, I, NP, IR, XI2, 
ETA2. ZETA2, SAI) 


Functional Description 

This subroutine constructs the integration matrix [S] described 
in Equation (3) a row at a time. 



SAIJ = 




for force calc, 
for moment calc. 



[ -cose- sinfj + cosy .• cosyj]/^ 


B.j = {cosa i [(y i -n i ^ 1 ^)(c0SYj) + (Zj-S^^sin^] 

- COSB. (Xj-C-j^kOSYj 

/ n k - va 

-cosY i (x j .-c i u; )sinY j }/c i 

where cosa.., cose.., cosy^ are the direction cosines of the input axis and 

where x., y., z., y- are the coordinates and dihedral of the aerodynamic box 
J J J J 

and are the coordinates of the one edge of the input axis. 

SAIJ is of course zero on boxes that are not to be integrated. 
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Description of Argument List 


FLAGA 


FLAGF 


XII 
ETA1 
ZETA1 
CG, SG 
CTIL 
X, Y, Z 
DELA 
LIMI 

UMAX 

I 

NP 

IR 

XI2 

ETA2 

ZETA2 

SAI 


Axis input option flag 
= 0, axis endpoints are input 
= 1, direction cosines are input 
Force/moment flag 

= 0, C. ^ e ' is a force in direction of axis 
lei 

= 1, C.. ' ' is a moment about axis 

’ Axis endpoint coordinates, 

Cosine-, sine of box dihedral angles 

Constant used to nondimensionalize integrated data ( 

Coordinates of the pressure points (boxes) 

Box areas 

First-, last box numbers for the integration of the 
ACp values 

Number of LIMI sets input for one constraint 

Intermediate index 

Number of ACp values 

Row index of the matrix 

Second axis endpoint coordinates when FLAGF = 0, 
direction cosines when FLAGF = 1 

A row of the integration matrix S 


Calling Subroutine SAIJ 
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SUBR0UTINE WEYKNP, NC. NEM, NELIMS. NM0N, NAXIS, NMIN, NMAX, NS, 

NPIT. NP0T. W, PI, DCP. EPS, PHI, SAI, DCPTIL, 

DELCPB, C0L, CIE, DCM0D. EBMIN, EBMAX, EB, ELIM, 

SBB, SBI , S, SBMAT, PCI. W0RK) 

Functional Description 

This subroutine is the core of the correction factor method. All logic 
for the method is established here. This subroutine uses input to decide 
what is to be done and sets up the argument lists for and executes the 
calls to all required subroutines. The following flow charts document the 
logic flow of this subroutine. This subroutine sets up the logic for 
various forms of input data and various types of calculations. The input 
data ranges over geometry, pressures, downwashes, aerodynamic influence 
matrices, previously generated correction factors, integration matrix data, 
etc. This data can enter the program by cards, tapes or both. 

There are three basic computational branches; (1) correction factor 
calculation, (2) monitoring of data (integration of pressures into aerodynamic 
parameters) and (3) application of previously generated correction factors 
to pressure distributions. Within branch (1) there exists a choice of what 
type of correction factors to generate, premultiplier, postmuliplier and 
new postmultiplier. Also a choice as to the type of weighting to be used 
(i.e. the T) is available. The program also tests to see If limits are placed 
on the correction factors and if modes are used. The constraining power a 
is always input since a constraint is simply a = 1.0. 
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Description of Argument List 


NP 

NC 

NEM 

NELIMS 

NM0N 

NAXIS 

NMIN 

NMAX 

NS 

NPIT 

NP0T 

W 

DI 

DCP 

EPS 

PHI 

SAI 

DCPTIL 

DELCPB 


Number of ACp elements 

Number of constraints 

Number of correction factor modes 

Number of input cards for the 

EBMIN, EBMAX pairs 

Number of sets of monitoring data 

Number of axes for use in the inegration of the 

AC values 
p t 

= max (1, NELIMS) 

= max (NC, NMIN, 10) 

= max (NP+NC, NEM+NC) 


Data set number of the system input data set 

Data set number of the system output data set 

A column of the normalwash matrix 

A row of the inverse-D matrix , (A matrix) 

A column of the theoretical AC matrix 

p t 

Incremental weight factors array, e 
A column of the weight factor mode shape matrix, + 
Integration matrix row array , [S]. SAN = [S] for monitoring 
A column of pressures modified by weight matrix, AC 

p e 

A column of the unmodified lifting pressure coefficients, ACp 
(either input, ACp, or computed as [D]" 1 {W}) 


t 
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C0L 

CIE 

DCM0D 

EBMIN, EBMAX 
EB 

ELIM 
SBB 
SB I 
S 

SBMAT 

DCI 

W0RK 


Complex array for intermediate use 
Array of the experimental constraints, C 

1 

Array of the modified AC e values 

Minimum-, maximum values allowed for the e array to take 
e array (incremental weight factors) (e = 4>e ) 

Array of the modified e values 
A row of the S matrix 
A row of the [S] [ >/T* matrix 

A two-dimensional complex work array of dimension NC by NC 
S matrix of maximum dimension NC by NS 
Array cf the AC g values 

A two-dimensional complex array of dimension NP by 



Calling Subroutine 
Called Subroutines 


MAIN 

CEMN, DCPB, DCPT, DELC, EDBL, EPSJ, 
GINV, M0DF, PHIJ, RECD, SAIJ, SDBL, 
WSWA 
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Flow Charts (Calculation of S, <f>, W) 
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Flow Charts (Accounting for limits) 



/Call M0DF 
/ computes 
\aC. 


■'mod 


/tall GINV 
/ computes 
N. general i zedy 
\inverse 










SUBROUTINE WSWA(NPIT, NP0T, NUTL1 , NTAPW, K0DE, NP, NCjgL , 

NMAX. NMSYM, NMASYM, W) 

Functional Description 

Subroutine WSWA is called from WEYT only if the input flag FLAGW = 1. 

It reads and prints the mode number and symmetry flag identifying the mode, 
the range of boxes over which the input W value applies, and the normal wash, 
M for this range of boxes. This card input is repeated for all ranges as 
needed, but only the non-zero W values are required as input. The complete 
W matrix is assembled from the input and is saved on tape NTAPW in column 

order. 
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Description of Argument List 


NPIT 

NP0T 

NUTL1 

NTAPW 

K0DE 

NP 

NC0L 

NMAX 

NMSYM 

NMASYM 

W 


Data set number of the system input data set 
Data set number of the system output data set 
Utility (scratch) tape number 
Tape number containing columns of the W matrix 
= -1 

Number of row elements in the W matrix 
Number of columns in the W matrix 
Maximum number of columns in the W matrix 
Number of symmetric modes 
Number of antisymmetric modes 

Two-dimensional complex array containing the W matrix 


Calling Subroutine WEYT 
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SUBR0UTINE ZER0UT(W0RK, LENGTH, L00P, ITAPE) 


Functional Description 

This subroutine initializes a complex array W0RK of length LENGTH 
to zeroes. In addition to this, when the argument ITAPE f 0, the complex 
zeroes stored in W0RK are written on tape ITAPE as many times as specified 
by the argument L00P. 


Description of Argument List 

W0RK Complex array to be initialized to zeroes 

LENGTH Length of the complex array W0RK 

L00P Number of times the array W0RK is to be written on 

tape ITAPE (only if ITAPE f 0) 

ITAPE Tape number on which the array W0RK is saved 

L00P-times (if any) 


Calling Subroutines EPSJ, MAIN, MATM, SDBL 
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■ Experiment [13] . 



x/c 


Figure 4. - Comparison of classic theory and two variations of 
the present method with experimental data for the 
steady case. 
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□ Experiment [13] 


Present Method fM^) 
Classic Theory 



Ftgur* 5. - Comparison of classic theory and the present method (N.) 

(with and without second order Bernoulli correction) with 

experimental data. 
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Figure 6. - Corape r Ison of various theories with experimental data for the 
244 oscillatory control surface case. 



Figure 7. - Comparison of classic theory and the present method 
(MJ with experimental data for the oscillatory case. 
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Figure 8. - Comparison of two variations of the Present transonic 

method with data and classic theory for the oscillatory case. 
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Figure 9. 



- Comparison of experimental and theoretical lift coefficient 
for an airfoil with a deflected 25* chord flap. 




Figure 10 - Comparison of experimental and theoretical pitching (about c/4) 
and hinge moment (about 3 c/4) coefficient for an airfoil 
with a deflected 25% chord flap. 
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Figure 12. - Comparison of experimental and theoretical pitching moment 
(about c/4) and hinge moment (about 3c/4) coefficient for 
an airfoil with an oscillating control surface. 
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Experiment k f = 0 
Uncorrected Theory 
Corrected Theory - 2 constraints 
Corrected Theory - 5 constraints 



Postmulti plying Correction Factors 


Experiment k r = 0 

Uncorrected Theory 

Corrected Theory - 2 constraints 
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Premultiplying Correction Factors 



Experiment k r = 0 

Uncorrected Theory 

Corrected Theory - 2 constraints 
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Postmulti plying Correction Factors 


Experiment 
Uncorrected Theory 
Corrected Theory - 2 constraints 
Corrected Theory - 5 constraints 



Using Premultiplying Correction Factors 




Experiment 
Uncorrected Theory 
Corrected Theory - 2 constraints 
Corrected Theory - 5 constraints 



Figure 19. - Real Part of Experimental, Theoretical and Corrected Pressure Loadings for Angle of Attack 
Using Postmultiplying Correction Factors 
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Using Postmultiplying Correction Factors 


Experiment 

Uncorrected Theory 

Corrected Theory - 2 constraints 



Premultiplying Correction Factors 



Experiment 

Uncorrected Theory 

Corrected Theory - 2 constraints 



Figure 22. - Imaginary Part of Experimental, Theoretical and Corrected Pressure Loadings for Flap Rotation 
Using Premultiplying Correction Factors 
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Figure 25. - Comparison of experimental data with corrected and uncorrected theory for a swept wing 
operating at subsonic speeds and a - 4°. 



Figure 26. - Comparison of experimental data with corrected ^nd uncorrected theory for a swept wing operating 
at subsonic speeds at a = 4°. 
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Figure 28. - Comparison of experimental data with corrected and uncorrected theory for a swept wing operating 
at subsonic speeds at a = 10°. 
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Figure 29. - Comparison of corrected and uncorrected classic theory and 
Present Method (M ) with experimental data. 




F1gur« 30. - Comparison of experimental deti with corrected theory; the cerr#ctf«B» r 
being based on static and oscillatory data. 



Figure 31. - Effect of applying static correction factors to the 
Present Method (M ) for the oscillatory case. 
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Flgure 32. - Effect of applying static correction factors to the 
Present Method (M l) for the oscillatory case. 
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Figure 33. - Application of static correction factors to the 
oscillatory case. 
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Deduced from Experiment [l 3 ] 
Correction Factor (New Post. 
Theory 



Figure 34. -Comparison of theoretical, corrected, ana experimental downwash induced by a control surface 
deflection on a two-dimensional section. 
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Figure 35. - Comparison of the apparent camber distribution induced at sub- and transonic 
Mach Numbers with the geometric camber for a control surface deflection on a 
two-dimensional section. 



o 



Figure 36. - Comparison of data with premultlpller and postmultiplier 
(New) corrected theory, 
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Figure 37. t Application of correction factors obtained at subsonic 

Mach Number (M = 0.5) to the transonic case (M^ = 0.85). 













Figure 40. 


Effect of applying correction factors obtained at = 0.5 
to other Mach Numbers for the Present Method (M ) 

oo ' 
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□ Experiment [13] 



Figure 41. - Results of applying postmultlplylng correction factors 
(New) obtained at M =0.5 and premultiplying factors 
obtained at M = 0.875 to theory. 
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Figure 43. - Comparison of experimental data with corrected and uncorrected theory for an arrow wing 
in pitch (a = 4°) operating at supersonic speeds. 





Figure 45. - Comparison of experimental data with corrected and uncorrected 
tneory for an Arrow wing In pitch. 
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Camber changes Inferred by New Postmultl plying correction factors 
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